10 REM       REFRACTION
20 REM
30 DEFSNG A-Z: REM <<<<<
40 INPUT "APPARENT ZENITH DIST";Z0
50 INPUT "AIR TEMP (F) ";T0
60 INPUT "REFR INDEX   ";N0
70 INPUT "HEIGHT (FEET)";H
80 R0=1: REM  PLANET RADIUS
90 M=1: REM  PLANET MASS
100 W=28.97: REM  MOL WT OF AIR
105 P=3.14159: P2=1.5708
106 RD=57.2957
110 REM
120 REM  CONVERT TO CGS UNITS
130 Z0=Z0/RD
140 H=H*30.48
150 T0=(T0-32)/1.8+273.1
160 W=W*1.665E-24
170 R0=R0*6.378E+08
180 M=M*5.976E+27
190 REM
200 REM  CALC ATMOS QUANTITIES AT R0
210 G=6.67E-08*M/(R0*R0)
220 S=1.38E-16*T0/(G*W)
230 LA=-.4*W*G/1.38E-16
240 GOSUB 350
250 REM   PRINT RESULTS
260 ZA=Z0*RD
270 ZR=(Z0+RF)*RD
280 RF=RF*3600*RD
290 PRINT "REFRACTION IS ";RF;" ARC SEC"
300 PRINT "PATHLENGTH IS ";AM;" AIR MASSES"
310 PRINT "APPAR ZENITH DIST ";ZA
320 PRINT "REAL ZENITH DIST  ";ZR
330 END
340 REM
350 REM  REFRACTION SUBR
360 REM
370 N=N0: DH=S/200: IF H<0 THEN DH=-DH
380 H1=0
390 T=T1: GOSUB 1110: T1=T
400 IF (H1-H)/DH>=0 THEN 450
410 H1=H1+DH: GOSUB 1110
420 EX=-1-(T0/(S*LR))
430 N=1+(N-1)*((T/T1)^EX)
440 T1=T: GOTO 400
450 REM  INITIALIZE PARAMETERS
460 RF=0: L=0: LZ=0: R=R0+H
470 L1=LR: Z=Z0
480 REM
490 REM  LAYER THICKNESS
500 IF R>R0+H+5*S THEN DR=S/10
510 IF R<=R0+H+5*S THEN DR=S/20
520 IF R<=R0+H+2*S THEN DR=S/50
530 IF R<=R0+H+.2*S THEN DR=S/200
540 F=1-DR/R: GOSUB 1170 : ZG=P-F
550 IF Z<=P2 THEN IN=1
560 IF Z>P2 THEN IN=0
570 IF Z>ZG THEN IN=-1
580 REM  INDEX OF REFR FOR SHELL
590 H1=R-R0: GOSUB 1110
600 T1=T-L1*DR*IN
610 EX=-1-(T0/(S*L1))
620 N=1+(N-1)*((T/T1)^EX)
630 L1=LR
640 IF IN=-1 THEN 690
650 IF IN=0 THEN 840
660 IF IN=1 THEN 940
670 REM
680 REM
690 REM  INWARD RAY
700 T1=T-LR*DR
710 N1=1+(N-1)*((T1/T)^EX)
720 T2=T-LR*DR*2
730 N2=1+(N-1)*((T2/T)^EX)
740 F=SIN(Z)*R/(R-DR): GOSUB 1170 : A=P-F
750 L=L+(N1-1)*R*SIN(Z-A)/SIN(A)
760 IF SIN(A)>N2/N1 THEN 800
770 F=SIN(A)*N1/N2: GOSUB 1170 : Z=P-F
780 RF=RF+(Z-A)
790 GOTO 820
800 Z=P-A
810 RF=RF-(P-2*Z)
820 GOTO 1050
830 REM
840 REM  GRAZING RAY
850 T1=T-LR*DR
860 N1=1+(N-1)*((T1/T)^EX)
870 A=P-Z
880 L=L+(N1-1)*(-2)*R*COS(Z)
890 F=SIN(A)*N1/N: IF F>1 THEN F=1
900 GOSUB 1170 : Z=F
910 RF=RF+(Z-A)
920 GOTO 1050
930 REM
940 REM  OUTGOING RAY
950 T2=T+LR*DR
960 N2=1+(N-1)*((T2/T)^EX)
970 F=SIN(Z)*R/(R+DR): GOSUB 1170 : A=F
980 L=L+(N-1)*R*SIN(Z-A)/SIN(A)
990 IF R>=R0+H THEN LZ=LZ+(N-1)*DR
1000 IF SIN(A)>N2/N THEN 1020
1010 F=SIN(A)*N/N2: GOSUB 1170 : Z=F: GOTO 1030
1020 Z=P-A
1030 RF=RF+(Z-A)
1040 REM
1050 REM  END SUBROUTINE
1060 R=R+DR*IN
1070 IF R<=R0+8*S THEN 490
1080 AM=L/LZ
1090 RETURN
1100 REM
1110 REM  TEMPERATURE SUBR
1120 LR=.5*LA: T=T0+LR*H1
1130 IF H1<=2*S THEN 1150
1135 LR=-.16*LA
1140 T=T0+S*LA+LR*(H1-2*S)
1150 RETURN
1160 REM
1170 REM   ARC SINE FUNCTION
1180 IF F<1 THEN 1200
1190 F=P2: GOTO 1210
1200 F=ATN(F/SQR(1-F*F))
1210 RETURN
1220 REM
1230 REM ************************************
1240 REM   APPEARED IN ASTRONOMICAL COMPUTING
1250 REM   SKY & TELESCOPE - MARCH 1989 ISSUE
1260 REM ************************************