home *** CD-ROM | disk | FTP | other *** search
/ ftp.barnyard.co.uk / 2015.02.ftp.barnyard.co.uk.tar / ftp.barnyard.co.uk / cpm / walnut-creek-CDROM / SIMTEL / CPMUG / CPMUG049.ARK / ASTRO.FOR < prev    next >
Text File  |  1984-04-29  |  6KB  |  231 lines
  1.     subroutine jayday(ny,nm,nd,jnl)
  2.     integer ny,nm,nd,my,mm
  3.         real jnl
  4.     my=ny
  5.     mm=nm
  6.     if (mm .gt. 2 ) go to 12
  7.     mm=mm + 9
  8.     my=my -1
  9.     go to 13
  10. 12    mm=mm-3
  11. 13      jnl =15078.0 + aint(1461.0*float(my)/4) 
  12.         jnl = jnl + aint((153.0*float(mm)+ 2.0)/5.0) + float(nd)
  13.     return
  14.     end
  15.     subroutine calday(jnl,ny,nm,nd)
  16.         real jnl
  17.     if ( jnl .lt. 15384.0 ) go to 10
  18.     if (jnl .lt. 88067.0) go to 11
  19. 10    nm=0
  20.     nd=0
  21.     ny=0
  22.     goto 13
  23. 11    nseq=jnl-15078
  24.     ny=(4*nseq-1)/1461
  25.     nd=(4*nseq + 3 -1461*ny)/4
  26.     nm= (5*nd-3)/153
  27.     nd=(5*nd+2-153*nm)/5
  28.     if ( nm .ge. 10) goto 12
  29.     nm=nm +3
  30.     goto 13
  31. 12    nm=nm-9
  32.     ny=ny+1
  33. 13    return
  34.     end
  35.     function sidtim(jday,fday)
  36.     double precision const0,const1,dthdt,a
  37.         real jday
  38.     dthdt=1.00273790929
  39.     twopi=6.283185308
  40.     const0=.277987616
  41.     const1 = .27379092965e-2
  42.     a= const0 + const1*dble(jday -33282.0)
  43.     a =a + dthdt*dble(fday)
  44.     sidtim= amod(sngl(a),1.0)*twopi
  45.     return
  46.     end
  47.     subroutine eapos(x,y,z,t)
  48.     real w,p,ec
  49.     real *8 mdt,m,twopi,dtu
  50. C     
  51. C ------Co-ordinates of earth in Ecliptic reference frame of date
  52. C ------ only approximate ( i.e. longitude accuracy +/- 1 minute of Arc
  53. C ------ Reference Almanac for computers -1980
  54.     data mdt,twopi/6.283019465D+2,6.28318531D+0/
  55.     dtu = dble(t-15019.5)/3.6525D+4
  56.     m= dble(-.026601727) + dtu*mdt
  57.     m = dmod(m,twopi)
  58.     w = sngl(m)
  59.     w = w +.0334409*sin(w) + 3.4907E-4*sin(2.0*w)
  60.     ec =.01675104 -.418E-4*sngl(dtu)
  61.     p= 1.0 -ec*ec
  62.     om = 1.766636813 + .030005*sngl(dtu)
  63.     xw =p*cos(w)/(1.0 + ec*cos(w))
  64.     yw = p*sin(w)/(1.0 + ec*cos(w))
  65.     co = cos(om)
  66.     so = sin(om)
  67.     x = xw*co -yw*so
  68.     y = xw*so + yw*co
  69.     z = 0
  70.     return
  71.     end
  72.     function dot(v1,v2)
  73.     dimension v1(3),v2(3)
  74. C -----  Vector dot product
  75.     dot = 0.0
  76.     do 1 i =1,3
  77. 1    dot = dot + v1(i)*v2(i)
  78.     return
  79.     end
  80.     subroutine cross(v1,v2,v3)
  81. C  
  82. C ----- Cross product of two vectors
  83.     dimension v1(3),v2(3),v3(3)
  84.     v3(1) = v1(2)*v2(3) -v1(3)*v2(2)
  85.     v3(2) =v1(3)*v2(1) - v1(1)*v2(3)
  86.     v3(3) = v1(1)*v2(2) -v1(2)*v2(1)
  87.     return
  88.     end
  89.      function decml(ch,n)
  90. C    DECML .. function for BCD to floating point binary conversion
  91. C    Nicole Simon
  92. C 
  93. C    Revised for FORTRAN-80 by Anthony BERESFORD 30 DEC 1979
  94.     logical*1 ch(1),bcd(13)
  95.     data bcd/1h1,1h2,1h3,1h4,1h5,1h6,1h7,1h8,1h9,1h-,1h.,1he/
  96.     data bcd(13)/1hE/
  97. C   
  98.      L=lstrng(ch,n)
  99. 6    in=n
  100. 7    kexpon = 0
  101. 8    kpoint = 0
  102. C    convert bcd to fixed point number
  103. 10    decml= 0
  104. 20    fmlplr = 10.**l
  105. 25    last = n+l -1
  106. 30    do 110 k= in,last
  107. 50    do 90 intger= 1,9
  108. 60    IF (ch(k) .ne. bcd(intger))goto 90
  109.     fmlplr = fmlplr/10.
  110.     decml= decml + float(intger)*fmlplr
  111.     goto 110
  112. 90    continue
  113.     if (ch(k) .eq. bcd(10))goto 100
  114.     if (ch(k) .eq. bcd(11))goto 105
  115.     if(ch(k) .eq. bcd(12) .or. ch(k) .eq. bcd(13))goto 115
  116.     goto 101
  117. 100    fmlplr = -fmlplr
  118. 101    fmlplr = fmlplr/10.
  119.     goto 110
  120. 105    kpoint = k
  121. 110    continue
  122. C   
  123. 111    goto 120
  124. C        
  125. 115    kexpon = k
  126.     in = k + 1
  127.     il = last -kexpon
  128. 120    if ( kpoint .ne. 0 ) goto 140
  129.     if (kpoint .eq. 0 .and. kexpon .eq. 0)goto 210
  130.     kpoint =kexpon
  131. 140    dvisor = 10.**(last-kpoint+1)
  132.     decml =decml/dvisor
  133. C       
  134.     if ( kexpon .eq. 0)goto 210
  135. C    convert to floating point number,shifting for exponent
  136.     kpoint = kdecml(ch,in,il)
  137.     decml =decml*10.**kpoint
  138. C       
  139. 210    n = n + l
  140. 220    return
  141.     end
  142.     function idecml(ch,n)
  143.     logical*1 ch(1)
  144. C        IDECML  .... function for bcd to fixed point binary conversion
  145. C        Nicole Simon
  146.     data limit/6/
  147. C    limit = maximum no. of digits allowed for integers
  148. 5    l =lstrng(ch,n)
  149.     if (l .lt. limit )goto 10
  150.     idecml =kdecml(ch,n,limit)
  151.     goto 220
  152. 10    idecml = kdecml(ch,n,l)
  153. 220    n= n + l
  154.        return
  155.     end
  156.     function kdecml(ch,n,l)
  157.     logical*1 ch(1),bcd(10)
  158.     real mulplr
  159.     data bcd/1h1,1h2,1h3,1h4,1h5,1h6,1h7,1h8,1h9,1h-/
  160. C        
  161.     kdecml = 0
  162.     mulplr = 10.**l
  163.     last = n+ l -1
  164.     do 110 k =n,last
  165.      mulplr =mulplr/10.
  166.     do 90 intger = 1,9
  167.     if (ch(k) .ne. bcd(intger))goto 90
  168.     kdecml = kdecml + ifix(float(intger)*mulplr)
  169.     goto 110
  170. 90    continue
  171.     if (ch(k) .ne. bcd(10) ) goto 110
  172.     mulplr = - mulplr
  173. 110    continue
  174.     return
  175.     end
  176.     function lstrng(ch,n)
  177. C       LSTRNG  ... function to search for length of input string
  178. C       starting in column  N
  179.     logical *1 ch(1),blank
  180.     data last,label,blank/80,81,1h /
  181.     do 50 k = n,last
  182.     if ( ch(k) .ne. blank )goto 50
  183.     lstrng = k -n
  184.     if( lstrng .ne. 0) goto 70
  185.     n = n +1
  186. 50    continue
  187.     lstrng = label -n
  188. 70    return
  189.     end
  190.     subroutine eaxyz(x,y,z,t,day)
  191. C * * * to give Earths's co-ordinates in Cartesian form
  192. C * * * for the reference frame 1950.0
  193.     real*4 w,ec,lam,l1950,kep,ld1950
  194.     real*8 mdt,m,twopi,dtu
  195.     data mdt,twopi/6.28301945717D+2,6.2831853071795864D+0/
  196.     data t1950,prec/33282.0,2.435637E-4/
  197.     data obli0,obli1,ecrot/.4093197,-2.27111E-4,2.28105E-6/
  198.     dtu =(dble(t-15019.) +dble(day-0.5))/3.6525D+4
  199.     ec = .01675104 -.418E-4*sngl(dtu)
  200.     m= dble(-.026601727) + dtu*mdt
  201.     m=dmod(m,twopi)
  202.     w = sngl(m)
  203.     ecan = kep(w,ec)
  204.     se = sin(ecan)
  205.     ce = cos(ecan)
  206.     r= 1.0 -ec*ce
  207.     upsil = atan2(sqrt(1.-ec*ec)*se,ce-ec)
  208.     om = 1.766636813 + .030005*sngl(dtu)
  209.     dy = (t-t1950)/365.25
  210.     lam = om + upsil
  211.     l1950 = lam - prec*dy
  212.     beta =-ecrot*dy*sin(lam +.091367)
  213.     xe = r*cos(l1950)*cos(beta)
  214.     ye =r*sin(l1950)*cos(beta)
  215.     ze =r*sin(beta)
  216.     obli =obli0 +obli1*sngl(dtu)
  217.     x= xe
  218.     y =ye*cos(obli) -ze*sin(obli)
  219.     z= ye*sin(obli) +ze*cos(obli)
  220.     return
  221.     end
  222.     real function Kep(am,ecc)
  223.     data eps/1.0E-8/
  224.     eni = am + ecc*sin(am)
  225. 10    de =(eni -ecc*sin(eni)-am)/(1.0 -ecc*cos(eni))
  226.     eni = eni -de
  227.     if (abs(de) .gt. eps)goto 10
  228.       kep = eni
  229.     return
  230.     end
  231.