

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C Program PLANETF
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
        PROGRAM XPLANET
        BYTE COMENT(25),FILNAM(25),OPTION,
     1  RETRO(8),BLAH
        LOGICAL HLP
        DOUBLE PRECISION T,AN1(9),XJD,CONV,GMT,XL,YL,XP,YP,XN,
     1  DAY,PI2,RULER(12),WEEK(7),WEEKR(7),
     2  YBES,TT
        DIMENSION EC0(9),EC1(9),AN0(9),P0(9),P1(9),EC(9),
     1  TL(9),TB(9),TH0(9),TH1(9),XI(9),R(9),A(9),G(9),GE(9),
     2  X(9),Y(9),Z(9),XX(9),YY(9),ZZ(9),IH(4,3),H(4),
     3  ZOD1(12),P(120,14),ID(120,14),IM(120,14),IS(120,14),
     4  GLYPH(120,14),IASP(14,14)
     8  ,MONTH(12),MOND(12),XNODE(2,10),YNODE(2,10)
CCCCC8CCCC
C
CCCCCCCCC
        DATA MOND/31,0,31,30,31,30,31,31,30,31,30,31/
        DATA MONTH/'JA','FE','MR','AP','MY','JN','JL','AU',
     1  'SE','OC','NO','DE'/
        DATA WEEK/'Monday  ','Tuesday ','Wednsday',
     1  'Thursday','Friday  ','Saturday','Sunday  '/
        DATA WEEKR/'Moon)   ','Mars)   ','Mercury)','Jupiter)',
     1  'Venus)  ','Saturn) ','Sun)    '/
        DATA ZOD/'Aries   ','Taurus  ','Gemini  ','Cancer  ',
     1  'Leo     ','Virgo   ','Libra   ','Scorpio ',
     2  'Sagtarus','Caprcorn','Aquarius','Pisces  '/
        DATA ZOD1/'A','T','G','C','L','V','=','S','/','K','Q','P'/
CCCCCCCCC
C
C     PLUTO ELEMENTS FROM SHARAF (1964)
C
CCCCCCCCC
        DATA A/1.00000023,.387098599,.723331619,1.523688395,
     1  5.202802875,9.53884320,19.19097811,30.0706724,39.672599/
        DATA AN0/358.475833,102.279381,212.603222,319.529425,
     1  225.444651,175.758444,74.313628,41.269550,231.002308/
        DATA AN1/35999.049750D0,149472.515289D0,58517.803875D0,
     1  19139.858500D0,3034.906654D0,1222.116782D0,428.502578D0,
     2  218.466783D0,144.072477D0/
        DATA EC0/.01675104,.20561421,.00682069,.09331290,
     1  .04825382,.05606075,.04704433,.00853341,.24706226/
        DATA EC1/-.00004180,.00002046,-.00004774,.00009206,
     1  .0,.0,.0,.0,.0/
        DATA P0/101.220833,75.899697,130.163833,334.218203,
     1  11.907422,90.110981,169.048778,43.755611,221.592475/
        DATA P1/1.719175,1.555489,1.408036,1.840758,
     1  .0,.0,.0,.0,1.388888/
        DATA TH0/0.,47.145944,75.779647,48.786442,
     1  98.932822,112.347606,73.490250,130.678889,108.937165/
        DATA TH1/0.,1.185208,.899850,.770992,
     1  .0,.0,.509667,1.100972,1.358056/
        DATA XI/0.,7.002881,3.393630,1.850333,
     1  1.311614,2.494239,.7726658,1.779256,17.109816/
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C Code
C
CCCCCCCCC
        TAN(DUMMY)=SIN(DUMMY)/COS(DUMMY)
        ARG(P,Q)=P-Q*INT(P/Q)-Q*INT(SIGN(.5,P-Q*INT(P/Q))-.5)
        HLP=.FALSE.
        LINEP=0
        PI=ATAN(1.0)
        PI2=6.283185307179586476D0
        CONV=PI2/360.D0
        CONVS=CONV
        DO 130 I=1,120
        DO 130 J=1,14
130     P(I,J)=0.
        WRITE(5,10000)
10000   FORMAT('$Sample input (Gregorian date,time',
     1  'conversion to GMT,latitude,longitude,OPTION)'//
     2  ' 1899/12/31 23:59 +08:00 34N59 122W59) : ')
CCCCCCCCC
C
C     HOROSCOPE LOOP STARTS HERE
C
CCCCCCCCC
        READ(5,11000)NY,CALEND,NM,ND,NH,NMI,NCH,NCM,LATD,LATS ,
     1  LATM,LONGD,LONGS,LONGM,OPTION
11000   FORMAT(I4,A1,4(I2,1X),I3,2(1X,I2),A1,I2,1X,I3,A1,I2,1X,1A)
6011    CONTINUE
        WRITE(6,11005)
11005   FORMAT('$Do you care for a hardcopy ? ')
        READ(5, 11006) BLAH
11006   FORMAT(A1)
        IF (BLAH.NE.'Y'.AND.BLAH.NE.'y') GO TO 11035
        CONTINUE
        CALL OPEN(6,'PLANET.DAT')
        IOUT=6
        GO TO 11008
11035   IOUT=5
11008   LINEP=1
        WRITE(IOUT,11500)NY,CALEND,NM,ND,NH,NMI,NCH,NCM,LATD,
     1  LATS,LATM,LONGD,LONGS,LONGM
11500   FORMAT('1'//' ---Date--- -Time -Zone- -Lat- -Long- ',
     1  'Comments--------'/' ',I4,A1,I2,'/',I2,I3,':',I2,
     1  I4,':',I2,I3,A1,I2,I4,A1,I2)
144     GMT=((NH+NCH-12)*60+NMI+SIGN(NCM*1.,NCH*1.))/1440.
        XJD=(ND-32075+1461*(NY+4800+(NM-14)/12)/4+367*(NM-2-
     1  (NM-14)/12*12)/12-3*((NY+4900+(NM-14)/12)/100)/4)+GMT
        IF(CALEND.NE.'B')GO TO 150
        NY=1-NY
150     IF(CALEND.NE.'J'.AND.CALEND.NE.'B')GO TO 160
        XJD=(ND-32075+1461*(NY+4800+(NM-14)/12)/4
     1  +367*(NM-2-(NM-14)/12*12)/12-38)+GMT
160     T=(XJD-2415020.D0)/36525.D0
C
        JFINAL=101
CCCCCCCCC
C
C     Planetary Ephemerides comxzted below (Newcomb/Hill,1898)
C
CCCCCCC
170     DO 270 J=1,JFINAL
        TT=T+(J-1)/36525.D0
        TP=TT+18262.D0/36525.D0
        W=TT
        DO 174 I=1,9
        EC(I)=EC0(I)+EC1(I)*W
        N=AN1(I)*TT/360.D0
174     G(I)=(AN1(I)*TT-N*360.D0+DBLE(AN0(I)))*CONV
        G5=G(5)
        G6=G(6)
        G7=(220.169542+428.49311*TP)*CONVS
        G(5)=G(5)+(.6506*SIN(2*G6-2*G5+336.9*CONVS)
     1  +(3.9987-.002213*36525./4332.58*TP)
     2  *SIN(5*G6-2*G5+(67.15-8197.0/3600.*TP)*CONVS)
     3  +.5380*SIN(5*G6-3*G5+176.5*CONVS)
     4  +.4112*SIN(2*G6-G5+1.4*CONVS)
     4  +.0399278*36525./4332.58*TP*SIN(-G5+227.46*CONVS)
     5  +.2763*SIN(3*G6-2*G5+127.4*CONVS)
     6  +.2669*SIN(G6-G5+79.2*CONVS)  )*299.12837/3600.*CONVS
CCCCCCCCC
C
CCCCCCCCC
        G(6)=G(6)+(24.153*SIN(5*G6-2*G5+(247.11-2.277*TP)*CONVS)
     1  +5.679*SIN(4*G6-2*G5+277.39*CONVS)
     2  +3.505*SIN(2*G6-G5+181.43*CONVS)
     3  +.657765*36525./10759.20*TP*SIN(G6+238.0*CONVS)
     4  +.278*SIN(3*G6-G5+121.2*CONVS)
     5  +.266*SIN(2*G6-2*G5+157.0*CONVS)
     6  +.238*SIN(6*G6-2*G5-3*G7+6.9*CONVS)
     7  +.223*SIN(10*G6-4*G5+(133.6-14814.5/3600.*W)*CONVS)
     8  +.234*SIN(3*G7-G6+321.7*CONVS)  )*120.455/3600.*CONVS
CCCCCCCCC
C
C
CCCCCCCCC
        G(7)=G(7)+33.086*W*W/3600.*CONVS
        G(8)=G(8)-22.401/3600.*W*W*CONVS
        DO 177 I=1,9
        GE(I)=G(I)+EC(I)*SIN(G(I)+EC(I)*SIN(G(I)+EC(I)*SIN(
     1  G(I)+EC(I)*SIN(G(I)+EC(I)*SIN(G(I))))))
177     X(I)=2.*ATAN(SQRT((1.+EC(I))/(1.-EC(I)))*TAN(.5*GE(I)))
     1  +(P0(I)+P1(I)*W)*CONVS + (I/7-I/9)*W*5025.3/3600.*CONVS
        GT=(358.415+35998.928*W)*CONVS
        GJ=(225.209+3034.462*W+.332*SIN((134.4+38.5*W)*CONVS))
     1  *CONVS
        X(4)=X(4)+(25.384*COS(GJ-G(4)-48.9*CONVS)
     1  +52.490*SIN((47.48+19.771*W)*CONVS)-37.05-13.50*W
     2  +21.869*COS(2*GJ-G(4)-188.3*CONVS)
     3  +16.035*COS(2*GJ-2*G(4)-191.9*CONVS)
     4  +13.966*COS(-GT+2*G(4)-20.5*CONVS)
     5  +8.559*COS(-GT+G(4)-35.1*CONVS)  )/3600.*CONVS
CCCCCCCCC
C
C
CCCCCCCCC
        GN=(225.417+3034.904*W)*CONVS
        GG=(175.753+1222.113*W)*CONVS
        G1=(74.412+428.498*W)*CONVS
        GP=(74.320+428.498*W)*CONVS
        GPP=(41.339+218.467*W)*CONVS
        X(7)=X(7)+(142.938*SIN(GG-2*G1)+19.508*COS(GG-2*G1)
     1  +75.70*COS(3*G1-GG)-102.30*SIN(3*G1-GG)
     2  -48.623*SIN(GN-G1)-21.320*COS(GN-G1)
     3  -27.871*COS(GP-GPP)+19.869*SIN(GP-GPP)
     4  +28.793*COS(2*GP-2*GPP)+10.035*SIN(2*GP-2*GPP)
     5  +(18.37*SIN(3*GP-3*GPP)+8.91*COS(3*GP-3*GPP))*COS(G(7))
     6  +(8.35*SIN(3*GP-3*GPP)-16.44*COS(3*GP-3*GPP))*SIN(G(7))
     7  -18.585*COS(GG-G1)+12.603*SIN(GG-G1)
     8  +4.327*COS(3*GP-3*GPP)+14.280*SIN(3*GP-3*GPP)
     9  )/3600.*CONVS
CCCCCCCCC
C
C
CCCCCCCCC
        X(7)=X(7)+((112.317*W-1.551*W*W-.516*W*W*W)*SIN(G(7))
     1  -(68.339*W+9.721*W*W)*COS(G(7))+6.605*W*SIN(2*G(7))
     2  -(29.44-.410*W)*SIN((20.45-22.61*W)*CONVS) )/3600.*CONVS
CCCCCCCCC
C
C     JUP/SAT TERM & QUADRATURES OF LONG PERIOD UR`WNEP TERMS
C
CCCCCCCCC
        X(8)=X(8)+(  33.972*W*COS(G(8))
     1  +18.553*SIN((180.966+1004.034*W+.1403*W*W)*CONVS)
     2  +34.138*SIN((153.267+2816.296*W-.0573*W*W)*CONVS)
     3  +(26.50*W+3.92*W*W)*SIN(G(8))   )/3600.*CONVS
CCCCCCCCC
C
C
CCCCCCCCC
        DO 185 I=1,9
        TH=(TH0(I)+TH1(I)*W)*CONVS
        DO 180 K=1,2
CCCCCCCCC
C
C     MOTION OF JUP/SAT NODES FROM LEVERRIER
C
CCCCCCCCC
        TH2=TH+(I/5*3636.6-I/6*493.1-I/7*3143.5)/3600.*TP*CONVS
        R(I)=A(I)*(1.-EC(I)*COS(2.*ATAN(TAN((TH2+(K-1)*PI-(P0(I)
     1  +P1(I)*W)*CONVS)*.5)*SQRT((1.-EC(I))/(1.+EC(I))))))
        XNODE(K,I)=R(I)*COS(TH2+(K-1)*PI)
180     YNODE(K,I)=R(I)*SIN(TH2+(K-1)*PI)
        TL(I)=TH+ATAN2(SIN(X(I)-TH)*COS(XI(I)*CONVS),
     1  COS(X(I)-TH))+(I/5-I/7)*TP*5026.1/3600.*CONVS
        TB(I)=ASIN(SIN(XI(I)*CONVS)*SIN(X(I)-TH))
185     R(I)=A(I)*(1.-EC(I)*COS(GE(I)))
        TB(5)=TB(5)+4.37431*36525./4332.58*SIN(X(5)+23.62*CONVS)
     1  *TP/3600.*CONVS
        TB(6)=TB(6)+24.266*36525./10759.2*SIN(X(6)-13.05*CONVS)
     1  *TP/3600.*CONVS
CCCCCCCCC
C
C     TB(5) EFFECTS LESS THAN 35"T, TB(6) LESS THAN 84"T
C
CCCCCCCCC
        R(5)=R(5)*10.**(.0002303*COS(2*G6-2*G5+336.9*CONVS)
     1  +.0001679*COS(5*G6-3*G5+176.4*CONVS)
     2  +.0000125634*36525./4332.58*TP*COS(-G5+227.4*CONVS)  )
CCCCCCCCC
C
C       22",16",10"T,(7",5",3")
C
CCCCCCCCC
        R(6)=R(6)*10.**(.0007005*COS(4*G6-2*G5+277.3*CONVS)
     1  +.0003783*COS(G6-G5+79.8*CONVS)
     2  +.000083491*36525./10759.20*TP*COS(G6+58.0*CONVS)
     3  +.0002443*COS(2*G6-G5+176.0*CONVS)  )
CCCCCCCCC
C
C
CCCCCCCCC
        DO 190 I=1,9
        X(I)=R(I)*COS(TB(I))*COS(TL(I))
        Y(I)=R(I)*COS(TB(I))*SIN(TL(I))
190     Z(I)=R(I)*SIN(TB(I))
        DO 210 I=2,9
        DO 200 K=1,2
        XNODE(K,I)=XNODE(K,I)-X(1)
        YNODE(K,I)=YNODE(K,I)-Y(1)
200     P(K+17,I)=ARG(ATAN2(YNODE(K,I),XNODE(K,I))/6.2831853,1.)
        XX(I)=X(I)-X(1)
        YY(I)=Y(I)-Y(1)
210     ZZ(I)=Z(I)-Z(1)
        XX(1)=-X(1)
        YY(1)=-Y(1)
        ZZ(1)=-Z(1)
        XNUT=-17.2327/1296000.*SIN((259.18-1934.142*W)*CONVS)
        DO 230 I=1,9
        P(J+19,I)=ARG(ATAN2(YY(I),XX(I))/6.28318531+XNUT,1.)
        IF(J-1)230,220,230
220     DIST=SQRT(XX(I)*XX(I)+YY(I)*YY(I)+ZZ(I)*ZZ(I))
        P(1,I)=ASIN(ZZ(I)/DIST)/CONVS
        E=(23.4522944-.0130125*W+.002558*COS((259.18-1934.142
     1  *W)*CONVS))*CONVS
        P(2,I)=ASIN(COS(P(1,I)*CONVS)*SIN(P(20,I)*6.28318531)
     1  *SIN(E)+SIN(P(1,I)*CONVS)*COS(E))/CONVS
        P(6,I)=ARG(TL(I)/6.28318531,1.)
        P(3,I)=TB(I)/CONVS
        P(4,I)=ASIN(COS(TB(I))*SIN(TL(I))*SIN(E)
     1  +SIN(TB(I))*COS(E))/CONVS
230     CONTINUE
CCCCCCCCC
C
C     LUNAR EPHEMERIDES COMPUTED TO WITHIN 1' (BROWN/I.L.E.)
C
CCCCCCCCC
        DAY=XJD-2415020.D0+(J-1.D0)
        XL=.751206D0+DAY*.0366011014634D0
        YL=.776935D0+DAY*.0027379092649D0
        XP=.928693D0+DAY*.0003094557786D0
        YP=.781169D0+DAY*.0000001307457D0
        XN=.719954D0-DAY*.0001470942283D0
        AL=(XL-IDINT(XL))*PI2
        BL=(YL-IDINT(YL))*PI2
        AP=(XP-IDINT(XP))*PI2
        BP=(YP-IDINT(YP))*PI2
        AN=(XN-IDINT(XN))*PI2
        U=AL-AP
        V=BL-BP
        F=AL-AN
        D=AL-BL
        DL=22639*SIN(U)-4586*SIN(U-D-D)+2370*SIN(D+D)+769*SIN(U+U)
     1  -668*SIN(V)-412*SIN(F+F)-212*SIN(U+U-D-D)-206*SIN(U+V-D-D)
     2  +192*SIN(U+D+D)-165*SIN(V-D-D)+148*SIN(U-V)-125*SIN(D)
     3  -110*SIN(U+V)-55*SIN(F+F-D-D)-45*SIN(U+F+F)+40*SIN(U-F-F)
     4  -38*SIN(U-4*D)+36*SIN(3*U)-31*SIN(U+U-4*D)+28*SIN(U-V-D-D)
     5  -24*SIN(V+D+D)+19*Sig(U-D)+18*SIN(V+D)+15*SIN(U-V+D+D)
     6  +14*SIN(4*D)+14*SIN(U+U+D+D)-13*SIN(3*U-D-D)
        P(J+19,10)=ARG(XL-IDINT(XL)+DL/1296000.+XNUT,1.)
        P(J+19,11)=ARG(AN/6.28318531,1.)
        IF(J-1)270,260,270
260     P(1,10)=(18461*SIN(F)+1010*SIN(U+F)-1000*SIN(F-U)
     1  -624*SIN(F-2*D)-167*SIN(U+F-D-D)+199*SIN(F-U+D+D)
     2  +117*SIN(F+D+D)+62*SIN(U+U+F)-33*SIN(F-U-D-D)
     3  -32*SIN(F-U-U)-30*SIN(V+F-D-D)-16*SIN(U+U+F-D-D)
     4  +15*SIN(U+F+D+D)+12*SIN(F-V-D-D)+9*SIN(F-U-V+D+D))/3600.
        P(2,10)=ASIN(COS(P(1,10)*CONVS)*SIN(P(20,10)*6.28318531)
     1  *SIN(E)+SIN(P(1,10)*CONVS)*COS(E))/CONVS
270     CONTINUE
CCCCCCCCC
C
CCCCCCCCC
        E=(23.4522944-.0130125*T+.002558*DCOS((259.18-1934.142
     1  *T)*CONV))*CONVS
        DO 277 I=2,9
        RETRO(I-1)=' DIR '
        IF(ARG(P(21,I)-P(20,I),1.)-.5)277,277,275
275     RETRO(I-1)='RETRO'
277     CONTINUE
CCCCCCCCC
C
CCCCCCCCC
        XLW=(LONGD+LONGM/60.)/360.
        XLN=(LATD+LATM/60.)*CONVS
        IF(LATS.EQ.'N')GO TO 280
        XLN=-XLN
280     IF(LONGS.EQ.'W')GO TO 290
        XLW=-XLW
CCCCCCCCC
C
C     H(1)=GMT  H(2)=GST  H(3)=LMT  H(4)=LST
C
CCCCCCCCC
290     H(1)=ARG(SNGL(GMT)+.5,1.)
        W=100.00213590D0*T-IDINT(100.D0*T)
        H(2)=ARG(H(1)+.27691940+W,1.)
        H(3)=ARG(H(1)-XLW,1.)
        H(4)=ARG(H(2)-XLW,1.)
        DO 300 I=1,4
        W=24.*H(I)+.00013888888
        IH(I,1)=W
        W=60.*(W-IH(I,1))
        IH(I,2)=W
300     IH(I,3)=60.*(W-IH(I,2))
        WRITE(IOUT,12000),((IH(I,J),J=1,3),I=1,4)
12000   FORMAT(/' GMT=',I2,2(':',I2),'    GST=',I2,2(':',I2),
     1  '    LMT=',I2,2(':',I2),'    LST=',I2,2(':',I2))
        JD1=XJD
        XJD1=XJD-JD1
        I=JD1+INT(XJD1+1.5-XLW)-(JD1-1+INT(XJD1+1.5-XLW))/7*7
        YBES=1900.D0+(XJD-2415020.31351528D0)/365.24219878125D0
     1  -.00000107523D0*T*T
        IYBES=YBES
        YBES1=YBES-IYBES
        WRITE(IOUT,13000)WEEK(I),WEEKR(I),JD1,XJD1,IYBES,YBES1
13000   FORMAT(' Day of the week (from LMT) IS ',A9,
     1  ' (ruled by ',A8/' JD=',I7,F7.6,7X,' Besselian year = ',
     2  I5,F10.9)
CCCCCCCCC
C
C     HOUSE CUSP COMPUTATIONS (ACCURATE TO LAST DIGIT)
C
CCCCCCCCC
310   ITER=1+LATD/45
CCCCCCCCC
C
C     "ITER" IS THE NUMBER OF PLACIDEAN ITERATIONS
C@@@
