1       REM
2       REM     Handling program for SUN
3       REM
5       DIM FLAG(7)
10      DEF FNRAD(W)=.01745329252#*W
15      DEF FNDEG(W)=57.29577951#*W
20      PRINT : FLAG(0)=1
25      PRINT "The position of the Sun"
30      PRINT : FLAG(1)=0 : FLAG(7)=0
35      INPUT "Date         (D,M,Y)"; DY,MN,YR
40      INPUT "... and time (H,M,S)"; XD,XM,XS
45      SW1=-1 : GOSUB 1000 : DY=DY+X/24
50      GOSUB 3300 : PRINT
55      PRINT "Geocentric coordinates:"
60      PRINT "---------- ------------"
65      PRINT : FLAG(6)=0
70      SW1=1 : X=FNDEG(LSN) : GOSUB 1000
75      PRINT "True longitude   (D,M,S): "; S$; XD; XM; XS
80      PRINT "True distance       (AU): "; RSN
85      X=LSN : Y=0 : SW3=-1 : GOSUB 1800
90      X=FNDEG(P/15) : GOSUB 1000
95      PRINT "True RA          (H,M,S): "; S$; XD; XM; XS
100     X=FNDEG(Q) : GOSUB 1000
105     PRINT "True declination (D,M,S): "; S$; XD; XM; XS
110     PRINT
115     GOSUB 2700 : LSND=FNDEG(LSN)+DPSI-.00569
120     X=LSND : GOSUB 1000
125     PRINT "Apparent longitude   (D,M,S): "; S$; XD; XM; XS
130     X=FNRAD(LSND) : Y=0 : SW3=-1 : GOSUB 1800
135     X=FNDEG(P/15) : GOSUB 1000
140     PRINT "Apparent RA          (H,M,S): "; S$; XD; XM; XS
145     X=FNDEG(Q) : GOSUB 1000
150     PRINT "Apparent declination (D,M,S): "; S$; XD; XM; XS
155     PRINT
160     INPUT "Again (Y/N)"; AN$
165     IF AN$="N" THEN STOP
170     GOTO 30
997     REM
998     REM     Subroutine MINSEC
999     REM
1000    IF SW1=-1 THEN GOTO 1035
1005    SN=SGN(X) : XP=ABS(X) : XD=INT(XP)
1010    A=(XP-XD)*60 : XM=INT(A)
1015    XS=INT((A-XM)*600+.5)/10
1020    S$="+"
1025    IF SN=-1 THEN S$="-"
1030    RETURN
1035    SN=+1
1040    IF XD<0 OR XM<0 OR XS<0 THEN SN=-1
1045    XD1=ABS(XD) : XM1=ABS(XM) : XS1=ABS(XS)
1050    X=((((XS1/60)+XM1)/60)+XD1)*SN
1055    RETURN
1097    REM
1098    REM     Subroutine JULDAY
1099    REM
1100    IF FLAG(1)=1 THEN RETURN
1105    DEF FNITG(W)=INT(W)+FLAG(0)*SGN(W)*INT((SGN(W)-1)/2)
1110    MN1=MN : YR1=YR : FLAG(1)=1 : B=0
1115    IF YR1<0 THEN YR1=YR1+1
1120    IF MN<3 THEN MN1=MN+12 : YR1=YR1-1
1125    IF YR<1582 THEN GOTO 1145
1130    IF YR=1582 AND MN<10 THEN GOTO 1145
1135    IF YR=1582 AND MN=10 AND DY<15 THEN GOTO 1145
1140    A=INT(YR1/100) : B=2-A+INT(A/4)
1145    IF YR1<0 THEN GOTO 1155
1150    C=INT(365.25*YR1)-694025! : GOTO 1160
1155    C=FNITG((365.25*YR1)-.75)-694025!
1160    D=INT(30.6001*(MN1+1))
1165    DJD=B+C+D+DY-.5
1170    RETURN
1697    REM
1698    REM     Subroutine OBLIQ
1699    REM
1700    IF FLAG(5)=1 THEN RETURN
1705    FLAG(5)=1
1710    GOSUB 1100 : T=DJD/36525!
1715    C=(((-.00181*T)+.0059)*T+46.845)*T
1720    EPS=23.45229444#-(C/3600)
1725    EPSR=EPS*.01745329252#
1730    RETURN
1797    REM
1798    REM     Subroutine EQECL
1799    REM
1800    IF FLAG(7)=1 THEN GOTO 1830
1805    DEF FNASN(W)=ATN(W/(SQR(1.000001-W*W)+9.999999E-21))
1810    PI=3.1415926535# : TPI=2*PI : FLAG(7)=1
1815    IF FLAG(6)=0 THEN DEPS=0
1820    GOSUB 1700 : EPS1=FNRAD(EPS+DEPS)
1825    SEPS=SIN(EPS1) : CEPS=COS(EPS1)
1830    CY=COS(Y) : SY=SIN(Y)
1835    IF ABS(CY)<9.999999E-21 THEN CY=9.999999E-21
1840    TY=SY/CY : CX=COS(X) : SX=SIN(X)
1845    SQ=(SY*CEPS)-(CY*SEPS*SX*SW3)
1850    Q=FNASN(SQ) : A=(SX*CEPS)+(TY*SEPS*SW3)
1855    P=ATN(A/CX)
1860    IF CX<0 THEN P=P+PI
1865    IF P>TPI THEN P=P-TPI
1870    IF P<0 THEN P=P+TPI
1875    RETURN
2697    REM
2698    REM     Subroutine NUTAT
2699    REM
2700    IF FLAG(6)=1 THEN RETURN
2705    REM DEF FNRAD(W)=1.745329252E-2*W
2710    FLAG(6)=1 : FLAG(7)=0
2715    GOSUB 1100 : T=DJD/36525! : T2=T*T
2720    A=100.0021358#*T : B=360*(A-INT(A))
2725    LS=279.697+.000303*T2+B
2730    A=1336.855231#*T : B=360*(A-INT(A))
2735    LD=270.434-.001133*T2+B
2740    A=99.99736056000026#*T : B=360*(A-INT(A))
2745    MS=358.476-.00015*T2+B
2750    A=13255523.59#*T : B=360*(A-INT(A))
2755    MD=296.105+.009192*T2+B
2760    A=5.372616667#*T : B=360*(A-INT(A))
2765    NM=259.183+.002078*T2-B
2770    TLS=2*FNRAD(LS) : NM=FNRAD(NM)
2775    TNM=2*FNRAD(NM) : MS=FNRAD(MS)
2780    TLD=2*FNRAD(LD) : MD=FNRAD(MD)
2785    DPSI=(-17.2327-.01737*T)*SIN(NM)+(-1.2729-.00013*T)*SIN(TLS)+.2088*SIN(TNM)-.2037*SIN(TLD)+(.1261-.00031*T)*SIN(MS)+.0675*SIN(MD)-(.0497-.00012*T)*SIN(TLS+MS)-.0342*SIN(TLD-NM)-.0261*SIN(TLD+MD)
2786    DPSI=DPSI+.0214*SIN(TLS-MS)-.0149*SIN(TLS-TLD+MD)+.0124*SIN(TLS-NM)+.0114*SIN(TLD-MD)
2790    DEPS=(9.21+.00091*T)*COS(NM)+(.5522-.00029*T)*COS(TLS)-.0904*COS(TNM)+.0884*COS(TLD)+.0216*COS(TLS+MS)+.0183*COS(TLD-NM)+.0113*COS(TLD+MD)-.0093*COS(TLS-MS)-.0066*COS(TLS-NM)
2795    DPSI=DPSI/3600 : DEPS=DEPS/3600
2800    RETURN
3097    REM
3098    REM     Subroutine ANOMALY
3099    REM
3100    TPI=6.283185308#
3105    M=MA-TPI*INT(MA/TPI) : EA=M
3110    DLA=EA-(S*SIN(EA))-M
3115    IF ABS(DLA)<.000001 THEN GOTO 3130
3120    DLA=DLA/(1-(S*COS(EA)))
3125    EA=EA-DLA : GOTO 3110
3130    TNU2=SQR((1+S)/(1-S))*TAN(EA/2)
3135    NU=2*ATN(TNU2)
3140    RETURN
3297    REM
3298    REM     Subroutine SUN
3299    REM
3300    GOSUB 1100 : T=DJD/36525! : T2=T*T
3305    REM DEF FNRAD(W)=1.745329252E-2*W
3310    A=100.0021359#*T : B=360*(A-INT(A))
3315    LS=279.69668#+.0003025*T2+B
3320    A=99.99736042000039#*T : B=360*(A-INT(A))
3325    MS=358.47583#-(.00015+.0000033*T)*T2+B
3330    S=.016751-.0000418*T-1.26E-07*T2
3335    MA=FNRAD(MS) : GOSUB 3100
3340    A=62.55209472000015#*T : B=360*(A-INT(A))
3345    A1=FNRAD(153.23+B)
3350    A=125.1041894#*T : B=360*(A-INT(A))
3355    B1=FNRAD(216.57+B)
3360    A=91.56766028#*T : B=360*(A=INT(A))
3365    C1=FNRAD(312.69+B)
3370    A=1236.853095#*T : B=360*(A-INT(A))
3375    D1=FNRAD(350.74-.00144*T2+B)
3380    E1=FNRAD(231.19+20.2*T)
3385    A=183.1353208#*T : B=360*(A-INT(A))
3390    H1=FNRAD(353.4+B)
3395    DL=.00134*COS(A1)+.00154*COS(B1)+.002*COS(C1)+.00179*SIN(D1)+.00178*SIN(E1)
3400    DR=5.43E-06*SIN(A1)+1.575E-05*SIN(B1)+1.627E-05*SIN(C1)+3.076E-05*COS(D1)+9.27E-06*SIN(H1)
3405    LSN=NU+FNRAD(LS-MS+DL) : TPI=6.283185308000017#
3410    RSN=1.0000002#*(1-S*COS(EA))+DR
3415    IF LSN<0 THEN LSN=LSN+TPI : GOTO 3415
3420    IF LSN>TPI THEN LSN=LSN-TPI : GOTO 3420
3425    RETURN