hp41programs

Refraction3

Ellipsoidal Refraction - Astronomical Refraction(3) for the HP-41


Overview
 

********************LATEST UPDATE*********************
*                                                                                                                *
*  This new version can be used for a triaxial ellipsoid model of the Earth     *
*                                                                                                                *
**********************************************************
 
 

-"EREF" calculates the refraction by rigorous differential equations.
-The 1st version used the classical 4th-order Runge-Kutta method to solve the system:

   dx/ds = A     dA/ds = ( n/x - A n/s ) / n                        x , y , z = cartesian coordinates of a running point
   dy/ds = B     dB/ds = ( n/y - B n/s ) / n              where      n     = refractive index of the Earth atmosphere
   dz/ds = C     dC/ds = ( n/z - C n/s ) / n                            ds    = element of the light ray curve

-But since it's actually a system of 3 second order differential equations, we can use other formulae to solve the system with less calculations, namely:

    y" = f(x,y,y')   is solved by:

    k = h f( x , y , y' )
   k' = h f( x+h/2 , y+(h/2).y'+(h/8).k , y'+k/2 )               and           y(x+h) ~ y(x) + h [ y' + ( k + k' + k" ) / 6 ]
   k" = h f( x+h/2 , y+(h/2).y'+(h/8).k , y'+k'/2 )                               y'(x+h) ~ y'(x) + ( k + k' + k" + k"' ) / 6
  k"' = h f( x+h , y+h.y'+(h/2).k" , y'+k" )
 

-The Earth may be approximated by an ellipsoid of revolution with semi-axes  a = 6378.137 km and  b = 6356.752314 km
  or by a triaxial ellipsoid with semi-axes  a = 6378.172 km , b = 6378.102 km , c = 6356.752314 km  and longitude of the major axis = 14°92911 W

>>>>   In the following program, clear F03 to employ an ellipsoid of revolution or set F03 to employ a triaxial ellipsoid.

-In the 1st version of "EREF" the distance point-ellipsoid was calculated by the following formula ( cf reference [3] ):

    Dist ~ E [ 1 / 4 / ( r2 / a4 + z2 / b4 ) + E ( r2 / a6 + z2 / b6 ) / 8 /  ( r2 / a4 + z2 / b4 )3 ] 1/2

    Where     E = ( r2 / a2 + z2 / b2 - 1 )   and   r2 = x2 + y2

-But I've found a better formula:

    Dist ~  [ ( x2 / a2 + y2 / b2 + z2 / c2 )1/2 - 1 ] [ a x2 + b y2 + c z2 + d ( x2 + y2 ) z2 ] / ( x2 + y2 + z2 )      with  d = - 6 E-9   if the length are expressed in km

-The errors don't exceed 10 mm over the interval  [-10 km , +100 km ]  , much less in a smaller interval.
-In fact, the errors would be even smaller if the calculator worked with more than 10 significant digits.
 

-The refractive index n is computed by   n = 1 + m  exp [ c1 Alt + c2 Alt2 + ............... + c13 Alt13 ]  ( 1 + k x + k' y + k" z2 )     with:
 

    c1 = t / 1250 - 0.109671                  c4 = -1.553002 E-4           c7 = 1.012373 E-8             c10 = 2529916 E-20            c13 = -41167039 E-28
    c2 = -0.0026952 - t x 95 E-7           c5 = 1.137826 E-5            c8 = -1054348 E-16           c11 = -31539538 E-23
    c3 =  958131 E-9                             c6 = -4.532222 E-7           c9 = -3737867 E-19           c12 = 1805402 E-24

  and  k = -1.585 E-6    k' = 1.718 E-6    k" = 2.855 E-9
 

>>>>      c1 to c13 are stored in R07 to R19  and  k , k' , k"  in R78 , R79 , R80.
 

-The constant "m" is deduced from the value of t P F Lambda which are to be stored in R01 thru R04.
-The refractivity at an altitude 0 ( or at the altitude of the observer ) is given by Owens' formula:

   ( n0 - 1 ) 108 = [ 2371.34 + 683939.7 ( 130 - 1/L2 ) -1 + 4547.3 ( 38.9 - 1/L2 ) -1 ] DS
                        + [ 6487.31 + 58.058 / L2 - 0.71150 / L4 + 0.08851 / L6 ] DW

   where    DS = ( PS / T ) [ 1 + PS ( 57.90 10 -8 - 9.325 10 -4 / T + 0.25844 / T2 ) ]
     and      DW = ( F / T ) [ 1 + F ( 1 + 0.00037 F ) ( -0.00237321 + 2.23366 / T - 710.792 / T2 + 77514.1 / T3 ) ]

   with   PS = pressure of dry air in mbar                     P = PS + F                              L = wavelength in µm
            F  = pressure of water vapor in mbar             T = temperature in Kelvin
 

>>>>  You can also use your own constants in registers R07 to R19 and R78 to R80 ( set F00 in this case )
 
 
 

Data Registers:           •  R00 = Longitude ( ° ' " ) positive East    ( Registers R00 thru R06 are to be initialized before executing "EREF" )

                                      •  R01 = t  ( °C )           •  R04 = wave-length ( µm )
                                      •  R02 = P ( mbar )       •  R05 = Latitude ( ° ' " )
                                      •  R03 = F ( mbar )       •  R06 = Altitude of the observer ( in meters )

                                          R07 = c1  ....................  R19 = c13      R20 to R84: temp   R51 = refractive index  R55 = increment/2   R77 = altitude

Flags:  F00 & F01         CF 00 = The constants in R07 to R19 & R78-79-80 are those listed below
                                        SF 00 = Store your own constants in R07 to R19 & R78-79-80 before executing "EREF"

                                       CF 01 = t , P , F  are measured for an altitude = 0
                                       SF 01 = t , P , F  are measured at the observer's station

                                       CF 03 = Ellipsoid of revolution
                                       SF 03 = Triaxial ellipsoid

                                       CF 04 means  k , k' , k" have the values given above if CF 00 or your own values in R78-R79-R80 if SF 00
                                       SF 04 means  k = k' = k" = 0
Subroutines: /
 

-Lines 236-516-517 are three-byte GTOs
 
 
 

 01 LBL "EREF"
 02 DEG
 03 HR
 04 STO 32
 05 X<>Y
 06 HR
 07 STO 33
 08 1
 09 P-R
 10 RCL 05
 11 HR
 12 STO 45            
 13 X<>Y
 14 P-R
 15 STO 36
 16 RDN
 17 CHS
 18 RCL 00
 19 HR
 20 14.92911
 21 +
 22 STO 66
 23 X<>Y
 24 P-R
 25 STO 34
 26 RDN
 27 STO 35
 28 CLX
 29 RCL 66
 30 X<>Y
 31 P-R
 32 ST+ 35
 33 X<>Y
 34 ST- 34
 35 RCL 32
 36 COS
 37 ST* 34
 38 ST* 35
 39 ST* 36
 40 RCL 45
 41 1
 42 P-R
 43 RCL 66
 44 X<>Y
 45 P-R
 46 RCL 32
 47 SIN
 48 ST* T
 49 ST* Z
 50 *
 51 ST+ 34
 52 RDN
 53 ST+ 35
 54 X<>Y
 55 ST+ 36
 56 6356.752314
 57 STO 65
 58 6378.137
 59 STO 41
 60 .035
 61 FC? 03
 62 CLX
 63 ST- 41
 64 +
 65 STO 40
 66 ST- Y
 67 RCL 65
 68 +
 69 *
 70 RCL 40
 71 X^2
 72 /
 73 RCL 45
 74 SIN
 75 X^2
 76 *
 77 RCL 41
 78 RCL X
 79 RCL 40
 80 ST- Z
 81 +
 82 *
 83 RCL 40
 84 X^2
 85 /
 86 RCL 45
 87 COS
 88 RCL 66
 89 SIN
 90 *
 91 STO 72
 92 X^2
 93 *
 94 +
 95 1
 96 +
 97 SQRT
 98 RCL 40
 99 X<>Y
100 /
101 STO 61
102 STO 62
103 STO 63
104 RCL 06
105  E3
106 /
107 STO 77
108 +
109 RCL 45
110 COS
111 RCL 66
112 COS
113 *
114 ST* 61
115 *
116 STO 37
117 RCL 41
118 RCL 40
119 /
120 X^2
121 ST* 62
122 RCL 63
123 *
124 RCL 77
125 +
126 RCL 72
127 ST* 62
128 *
129 STO 38
130 RCL 63
131 RCL 65
132 RCL 40
133 /
134 X^2
135 ST* 63
136 *		 137 RCL 77
138 +
139 RCL 45
140 SIN
141 ST* 63
142 *
143 STO 39
144 RCL 34
145 STO 42
146 RCL 35
147 STO 43
148 RCL 36            
149 STO 44
150 RCL 01
151 273.15
152 +
153 1/X
154 ENTER
155 ENTER
156 STO 53
157 77514.1
158 *
159 710.792
160 -
161 *
162 2.23366
163 +
164 *
165 .00237321
166 -
167 RCL 03
168 ST* Y
169 37 E-5
170 *
171 1
172 +
173 *
174 1
175 +
176 RCL 03
177 *
178 *
179 STO 51
180 RCL 04
181 X^2
182 1/X
183 ENTER
184 ENTER
185 STO 52
186 .08851
187 *
188 .7115
189 -
190 *
191 58.058
192 +
193 *
194 6487.31
195 +
196 ST* 51
197 RCL 53
198 0.25844
199 *
200 9325 E-7
201 -
202 RCL 53
203 *
204 579 E-9
205 +
206 RCL 02
207 RCL 03
208 -
209 ST* 53
210 *
211 1
212 +
213 RCL 53
214 *
215 4547.3
216 38.9
217 RCL 52
218 -
219 /
220 130
221 RCL 52
222 -
223 683939.7
224 X<>Y
225 /
226 +
227 2371.34
228 +
229 *
230 RCL 51
231 +
232  E8
233 /
234 STO 76
235 FS? 00
236 GTO 00
237 CLX
238 FC? 04
239 -1585 E-9
240 STO 78
241 FC? 04
242 1718 E-9
243 STO 79
244 FC? 04
245 2.855 E-9
246 STO 80
247 RCL 01
248 1250
249 /
250 .109671
251 -
252 STO 07
253 .0026952
254 CHS
255 RCL 01
256 95 E-7
257 *
258 -
259 STO 08
260 958131 E-9
261 STO 09
262 1.553002 E-4
263 CHS
264 STO 10
265 1.137826 E-5
266 STO 11
267 4.532222 E-7
268 CHS
269 STO 12
270 1.012373 E-8
271 STO 13
272 1054348 E-16		 273 CHS
274 STO 14
275 3737867 E-19
276 CHS
277 STO 15
278 2529916 E-20
279 STO 16
280 31539538 E-23
281 CHS
282 STO 17
283 1805402 E-24
284 STO 18
285 41167039 E-28
286 CHS
287 STO 19
288 LBL 00
289 31.019
290 STO 75
291 19
292 13
293 LBL 07
294 RCL IND Y
295 RCL Y
296 *
297 STO IND 75
298 CLX
299 SIGN
300 ST- Z
301 -
302 DSE 75
303 GTO 07
304 RCL 76
305 FC? 01
306 GTO 00
307 RCL 37
308 STO 61
309 RCL 38
310 STO 62
311 RCL 39
312 STO 63
313 RCL 77
314 XEQ 04
315 LASTX
316 LBL 00
317 RCL 61
318 RCL 78
319 *
320 RCL 62
321 RCL 79
322 *
323 +
324 RCL 63
325 X^2
326 RCL 80
327 *
328 +
329 1
330 +
331 *
332 RCL 76
333 /
334 ST/ 76
335 2
336 STO 82
337 6 E-9
338 CHS
339 STO 81
340 41
341 STO 84
342 LBL 01
343 RCL 77
344 RCL 84
345 /
346 E^X
347 STO 55
348 RCL 37
349 STO 45
350 RCL 38
351 STO 46
352 RCL 39
353 STO 47
354 RCL 42
355 STO 59
356 RCL 43
357 STO 60
358 RCL 44
359 STO 61
360 XEQ 02
361 STO 64
362 STO 70
363 STO 73
364 RCL 82
365 /
366 RCL 61
367 +
368 RCL 55
369 *
370 ST+ 47
371 RCL 75
372 STO 63
373 STO 69
374 STO 72
375 RCL 82
376 /
377 RCL 60
378 +
379 RCL 55
380 *
381 ST+ 46
382 RCL 74
383 STO 62
384 STO 68
385 STO 71
386 RCL 82
387 /
388 RCL 59
389 +
390 RCL 55
391 *
392 ST+ 45
393 RCL 62
394 ST+ 59
395 RCL 63
396 ST+ 60
397 RCL 64
398 ST+ 61
399 XEQ 02
400 ST+ 70
401 ST+ 73
402 ST+ 73
403 RCL 44
404 +
405 STO 61
406 RCL 75
407 ST+ 69
408 ST+ 72		 409 ST+ 72
410 RCL 43
411 +
412 STO 60
413 RCL 74
414 ST+ 68
415 ST+ 71
416 ST+ 71
417 RCL 42
418 +
419 STO 59
420 XEQ 02
421 STO 64            
422 ST+ 70
423 ST+ 73
424 ST+ 73
425 RCL 44
426 +
427 RCL 55
428 ST+ X
429 STO 83
430 *
431 RCL 39
432 +
433 STO 47
434 RCL 75
435 STO 63
436 ST+ 69
437 ST+ 72
438 ST+ 72
439 RCL 43
440 +
441 RCL 83
442 *
443 RCL 38
444 +
445 STO 46
446 RCL 74
447 STO 62
448 ST+ 68
449 ST+ 71
450 ST+ 71
451 RCL 42
452 +
453 RCL 83
454 *
455 RCL 37
456 +
457 STO 45
458 RCL 42
459 RCL 62
460 ST+ X
461 +
462 STO 59
463 RCL 43
464 RCL 63
465 ST+ X
466 +
467 STO 60
468 RCL 44
469 RCL 64
470 ST+ X
471 +
472 STO 61
473 XEQ 02
474 ST+ 73
475 RCL 75
476 ST+ 72
477 RCL 74
478 ST+ 71
479 PI
480 INT
481 ST/ 68
482 ST/ 69
483 ST/ 70
484 ST/ 71
485 ST/ 72
486 ST/ 73
487 RCL 42
488 RCL 68
489 +
490 RCL 83
491 *
492 ST+ 37
493 RCL 43
494 RCL 69
495 +
496 RCL 83
497 *
498 ST+ 38
499 RCL 44
500 RCL 70
501 +
502 RCL 83
503 *
504 ST+ 39
505 RCL 71
506 ST+ 42
507 RCL 72
508 ST+ 43
509 RCL 73
510 ST+ 44
511 VIEW 77
512 CLX
513 SIGN
514 RCL 51
515 X#Y?
516 GTO 01
517 GTO 12
518 LBL 02
519 RCL 45
520 RCL 40
521 /
522 X^2
523 RCL 46
524 RCL 41
525 /
526 X^2
527 +
528 RCL 47
529 RCL 65
530 /
531 X^2
532 +
533 SQRT
534 ENTER
535 STO 53
536 SIGN
537 -
538 STO 52
539 RCL 40
540 RCL 45
541 X^2
542 STO 49
543 *
544 RCL 41		 545 RCL 46
546 X^2
547 ST+ 49
548 *
549 +
550 RCL 65
551 RCL 47
552 X^2
553 STO 57
554 *
555 +
556 RCL 49            
557 ST+ 57
558 RCL 47
559 X^2
560 *
561 RCL 81
562 *
563 +
564 STO 58
565 RCL 57
566 /
567 STO 74
568 *
569 XEQ 04
570 STO 51
571 RCL 47
572 RCL 58
573 RCL 53
574 /
575 STO 58
576 *
577 RCL 65
578 X^2
579 /
580 RCL 49
581 RCL 81
582 *
583 RCL 65
584 +
585 RCL 74
586 -
587 RCL 47
588 ST+ X
589 *
590 RCL 52
591 *
592 +
593 STO 50
594 RCL 45
595 RCL 58
596 *
597 RCL 40
598 X^2
599 /
600 RCL 40
601 RCL 47
602 X^2
603 RCL 81
604 *
605 STO 75
606 +
607 RCL 74
608 -
609 RCL 45
610 ST+ X
611 *
612 RCL 52
613 *
614 +
615 STO 48
616 RCL 46
617 RCL 58
618 *
619 RCL 41
620 X^2
621 /
622 RCL 41
623 RCL 75
624 +
625 RCL 74
626 -
627 RCL 46
628 ST+ X
629 *
630 RCL 52
631 *
632 +
633 RCL 54
634 RCL 57
635 /
636 ST* 48
637 ST* 50
638 *
639 STO 49
640 FS? 04
641 GTO 03
642 RCL 45
643 RCL 78
644 *
645 RCL 46
646 RCL 79
647 *
648 +
649 RCL 47
650 X^2
651 RCL 80
652 *
653 +
654 PI
655 SIGN
656 +
657 ST* 48
658 ST* 49
659 ST* 50
660 RCL 67
661 *
662 PI
663 SIGN
664 +
665 STO 51
666 RCL 78
667 RCL 67
668 *
669 ST+ 48
670 RCL 79
671 LASTX
672 *
673 ST+ 49
674 RCL 80
675 ST+ X
676 LASTX
677 *
678 RCL 47
679 *
680 ST+ 50
681 LBL 03		 682 RCL 48
683 STO Y
684 RCL 59
685 *
686 RCL 49
687 RCL 60
688 *
689 +
690 RCL 50
691 RCL 61
692 *
693 +
694 STO 52            
695 RCL 59
696 *
697 -
698 STO 74
699 RCL 49
700 RCL 52
701 RCL 60
702 *
703 -
704 STO 75
705 RCL 50
706 RCL 52
707 RCL 61
708 *
709 -
710 RCL 55
711 RCL 51
712 /
713 ST* 74
714 ST* 75
715 *
716 RTN
717 LBL 04
718 ENTER
719 ENTER
720 STO 77
721 RCL 31
722 *
723 RCL 30
724 +
725 *
726 RCL 29
727 +
728 *
729 RCL 28
730 +
731 *
732 RCL 27
733 +
734 *
735 RCL 26
736 +
737 *
738 RCL 25
739 +
740 *
741 RCL 24
742 +
743 *
744 RCL 23
745 +
746 *
747 RCL 22
748 +
749 *
750 RCL 21
751 +
752 *
753 RCL 20
754 +
755 *
756 RCL 07
757 +
758 STO 54
759 CLX
760 RCL 19
761 *
762 RCL 18
763 +
764 *
765 RCL 17
766 +
767 *
768 RCL 16
769 +
770 *
771 RCL 15
772 +
773 *
774 RCL 14
775 +
776 *
777 RCL 13
778 +
779 *
780 RCL 12
781 +
782 *
783 RCL 11
784 +
785 *
786 RCL 10
787 +
788 *
789 RCL 09
790 +
791 *
792 RCL 08
793 +
794 *
795 RCL 07
796 +
797 *
798 E^X
799 RCL 76
800 *
801 ST* 54
802 SIGN
803 LASTX
804 STO 67
805 +
806 RTN
807 LBL 12
808 RCL 35
809 RCL 44
810 *
811 RCL 36
812 RCL 43
813 *
814 -
815 X^2
816 RCL 36
817 RCL 42
818 *		 819 RCL 34
820 RCL 44
821 *
822 -
823 X^2
824 +
825 RCL 34
826 RCL 43
827 *
828 RCL 35
829 RCL 42            
830 *
831 -
832 X^2
833 +
834 RCL 42
835 X^2
836 RCL 43
837 X^2
838 RCL 44
839 X^2
840 +
841 +
842 /
843 ABS
844 SQRT
845 ASIN
846 STO 52
847 RCL 66
848 RCL 33
849 SIN
850 STO 50
851 RCL 05
852 HR
853 STO 49
854 SIN
855 *
856 P-R
857 RCL 66
858 RCL 33
859 COS
860 P-R
861 RDN
862 -
863 RCL 42
864 *
865 X<> Z
866 +
867 RCL 43
868 *
869 +
870 RCL 50
871 RCL 49
872 COS
873 STO 53
874 *
875 RCL 44
876 *
877 -
878 STO 47
879 RCL 66
880 RCL 53
881 P-R
882 RCL 42
883 *
884 X<>Y
885 RCL 43
886 *
887 +
888 RCL 49
889 SIN
890 RCL 44
891 *
892 +
893 STO 46
894 RCL 66
895 RCL 49
896 RCL 33
897 COS
898 P-R
899 RCL 44
900 *
901 STO 45
902 RDN
903 P-R
904 RCL 66
905 RCL 33
906 SIN
907 P-R
908 R^
909 -
910 RCL 43
911 *
912 ST+ 45
913 RDN
914 +
915 RCL 42
916 *
917 ST- 45
918 RCL 32
919 1
920 P-R
921 RCL 46
922 *
923 CHS
924 X<>Y
925 RCL 45
926 *
927 +
928 RCL 45
929 X^2
930 RCL 46
931 X^2
932 +
933 X#0?
934 /
935 ASIN
936 STO 53
937 RCL 47
938 ENTER
939 X^2
940 RCL 45
941 X^2
942 +
943 SQRT
944 X#0?
945 /
946 ASIN
947 RCL 52
948 X<>Y
949 RCL 53
950 3600
951 ST* T
952 ST* Z
953 *
954 CLD
955 END

 
      ( 1714 bytes / SIZE 085 )
 
 

      STACK        INPUTS      OUTPUTS
           Z              /        RTotal
           Y         Az ( ° ' " )        RAzim
           X         h0 ( ° ' " )        RVertical

    All the refractions are expressed in arcseconds

Example:       t = 10°C  P = 1010 mbar  F = 6 mbar  Lambda = 0.577 µm  Latitude = 33°21'22" N  Altitude = 1706 m    Longitude = 116°51'50" W

    10    STO 01        0.577    STO 04
  1010  STO 02      33.2122  STO 05         -116.5150  STO 00
     6     STO 03        1706     STO 06

  1°)     CF 00  CF 01  CF 03  SF 04    Az = 12°41'   h0 = 1°23'45"

       12.41    ENTER^
      1.2345   XEQ "EREF"  >>>>  RV = 1076"3761                        ---Execution time = 17s---   with V41 in turbo mode.
                                           RDN   RAz = 0"0249
                                          X<>Y  RT = 1076"3762

  2°)     CF 00  CF 01  CF 03  CF 04    Az = 12°41'   h0 = 1°23'45"

       12.41    ENTER^
      1.2345       R/S       >>>>  RV = 1079"0328                   ---Execution time = 19s---   with V41 in turbo mode.
                                     RDN   RAz = -0"0263
                                    X<>Y  RT = 1079"0328
 

  3°)     CF 00  SF 01  CF 03  SF 04    Az = 12°41'   h0 = 1°23'45"

       12.41    ENTER^
      1.2345      R/S       >>>>  RV = 1295"9029
                                    RDN   RAz = 0"0301
                                   X<>Y  RT = 1295"9029

  4°)     CF 00  SF 01  CF 03  CF 04    Az = 12°41'   h0 = 1°23'45"

       12.41    ENTER^
      1.2345      R/S       >>>>  RV = 1299"1188
                                    RDN   RAz = -0"0314
                                   X<>Y  RT = 1299"1188

-If you set F03 for a triaxial ellipsoid, the differences are very small in the above example:
-For instance,

                    CF 00  CF 01  SF 03  CF 04

       12.41    ENTER^
      1.2345      R/S        >>>>  RV = 1079"0338  instead of  1079"0328
                                     RDN   RAz = -0"0262 ( almost unchanged )
                                    X<>Y  RT = 1079"0338  instead of  1079"0328

-With  Az = 80°  and  h0 = 0

           80  ENTER^
            0      R/S           >>>>  RV  = 1725"4899  instead of  1725"4763
                                     RDN   RAz = -0"2184    instead of  -0"2184
                                    X<>Y  RT = 1725"4900  instead of  1725"4763

-So, the difference is just about than 0"014
 

Notes:

-The azimuths are reckoned clockwise from North and the longitudes are positive East.

-With free42  execution time becomes less than 2 seconds !
-But with a real HP41, it may require about 6 hours...

-The increment is  2 Exp(Alt/41)  km   ( lines 343 to 347 )
-So,  2 km if Alt = 0  and about  16 km  if  Altitude = 100 km.

-Longitudes are positive East, latitudes are positive North.

-If F04 is set ( k = k' = k" = 0 ) and if F03 is clear, the longitude doesn't matter.
-In this case store, for instance,  0 in R00

-Horizontal refraction RAzim is ( almost ) always very small.

-The HP41 displays the successive altitudes.
-"EREF" stops when the refractive index = 1 ( line 515 )
 
 
 
 

References:

[1]  Exact Theory of Astronomical Refraction for the Real Atmosphere - A. Yu. Yatsenko - Astronomy & Astrophysics 579-586 ( 1995 )
[2]  https://ccmc.gsfc.nasa.gov/modelweb/models/msis_vitmo.php
[3]  Approximate formulas for the Point-to-Ellipse and for the Point-to-Ellipsoid Distance Problem - Alexei Uteshev & Marina Yashina