Announcement

Collapse
No announcement yet.

Four Dimensional Trig makes my brain hurt

Collapse
X
 
  • Filter
  • Time
  • Show
Clear All
new posts

  • #21
    Eric,
    try this version and see if it's more what you expect.
    All I've done is derive the angle that the satellite has moved from the equation for latitude instead of the equation for longitude.

    It now doesn't handle the quadrants as well but you said you were ok with that.
    Previously, numbers could not go out of range but now they can so you might want to look into range checking, maybe along the lines of what Mark posted.

    This gives a heading based on latitude.

    Your house is at 45N so, looking at the program output, the heading when passing over your house is 68.96degrees. You'll need to adjust for North or South.
    Does this look right?

    Code:
    'PBCC6 program
     #COMPILE EXE
     #BREAK ON
    
    
     GLOBAL Pi AS EXT
    
     FUNCTION PBMAIN
    
     LOCAL Angle,Lat,Lon,Cos516, Sin516, DegToRad, RadToDeg, Degrees, IncrementalLat, IncrementalLon, IncrementalAngle, Increment, Heading AS EXT
    
     Pi = 3.14159265358979323##
    
     Increment = 0.0000001##  'a very small increment to apply to the orbital position to determine which direction the satellite is moving
    
     'Factors to convert degrees to/from radians
     DegToRad = pi/180
     RadToDeg = 180/pi
    
    
     'The inclination of the orbit is a constant of 51.6 and it's useful to calculate these 2 values in advance as they don't change
     Cos516 = COS(DegToRad * 51.6)
     Sin516 = SIN(DegToRad * 51.6)
    
    
     'Angle is the distance around the orbit that the satellite has moved.
     'Zero degrees is the point over the equator, heading North.
     'Heading is the Angle from North, so North would be zero degrees heading, East is 90 degrees heading.
    
    
     PRINT "Degrees","Latitude","Heading"
    
    
     FOR Degrees = -52 TO 52  STEP 0.1##   'NOTE this is now the given LATITUDE, not the angle that the satellite has moved around its orbit
         'convert degrees to Radians
         Lat = Degrees * DegToRad
    
         'work out an angle that the satellite has moved through which corresponds to this latitude
         Angle = ArcSin(SIN(lat)/sin516)   + 0.000001
    
    
         'work out the latitude and longitude at this point in the orbit
         Lon = ATN2(COS(Angle), SIN(Angle)*Cos516)
         Lat = ArcSin(SIN(Angle)*Sin516)      'this latitude should be the same as the given latitude so it's not necessary to calculate it
    
         'To find the heading, move the satellite forward a small increment and recalculate the Latitude and Longitude
         IncrementalAngle = Angle + increment
    
         IncrementalLon = ATN2(COS(IncrementalAngle), SIN(IncrementalAngle)*Cos516)
         IncrementalLat = ArcSin(SIN(IncrementalAngle)*Sin516)
    
         'The heading is given by the direction the satellite has moved over the small increment in angle.
         Heading =  ATN2(IncrementalLat - Lat, IncrementalLon - Lon)
    
    
         PRINT FORMAT$(Degrees,"0.00"), FORMAT$((Lat*RadToDeg),"0.000"), FORMAT$(Heading*RadToDeg,"0.000")
    
     NEXT
    
    
     PRINT "done"
     WAITKEY$
    
    
     END FUNCTION
    
    
    
     FUNCTION ArcSin(Value AS EXT) AS EXT
    
     FUNCTION = ATN(Value / SQR(1 - Value * Value))
    
     END FUNCTION
    
    
    
    
     FUNCTION ArcCos(Value AS EXT) AS EXT
    
     FUNCTION = pi / 2 - ATN(Value / SQR(1 - Value * Value))
    
    
     END FUNCTION
    
    
    
     FUNCTION ATN2 ( X AS EXT, Y AS EXT) AS EXT
     #REGISTER NONE
    
     '!fld y                  ;Get y onto the FPU stack (for parameters passed ByVal)
     !mov eax,y              ;for ByRef (a little quicker)get the pointer to Y
     !fld tbyte ptr [eax]    ;load the  value pointed to by eax
    
    
     '!fld x                 ;Get x onto the FPU stack (for parameters passed ByVal)
     !mov eax,x              ;for ByRef (a little quicker)get the pointer to x
     !fld tbyte ptr [eax]    ;load the  value pointed to by eax
    
     !fpatan                 ;Get the ATAN(y/x), Put result In ST(1) Then pop the stack.
                             'i.e. result is now on top of stack in st(0)
    
     !ftst                   ;test ST(0) (that's the top of stack) against zero
     !fstsw ax               ;get the result into AX
     !sahf                   ;get the FPU flags into the CPU flags
     !jae short skip         ;if above or equal (i.e.value is not -ve) then skip the Add 2*pi code
    
     !fldpi                  ;get pi
     !fadd st(1),st(0)       ;add pi to result
     !faddp                  ;and again, for 2pi then pop the now unneeded pi from st(0)
                             'this leaves the result +2pi on the stack at st(0) if it was -ve
      skip:
    
     !fstp function          ;store the result directly in the function return value and pop the stack
    
    
     END FUNCTION

    Comment


    • #22
      Wow you're fast! I'll take a close look this weekend.

      And thanks again!
      "Not my circus, not my monkeys."

      Comment

      Working...
      X