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