module astro
!-------------------------------------------------------------
!  Sekalaisia pallotahtitieteen rutiineita
!  - juliaaninen paivamaara JD (standardi ja modifioitu) 
!  - paivan numero (oma JD:n kokonaislukuversio) ja vastaava
!    Gregoriaanisen kalenterin pvm
!  - tahtiaika
!  - Auringon ja Kuun apparentit ekvatoriaaliset koordinatit
!-------------------------------------------------------------

implicit none
private

double precision, parameter :: two_pi = 6.283185307
double precision, dimension(0:34) :: C
integer :: zone = 0 

public :: sidt0, sidt, & 
          julian, julian0, day_number, gregorian, &
          inits, sun, moon

contains

double precision function dmod(x, p)
  double precision, intent(in) :: x, p
  dmod = x-floor(x/p)*p
end function

!-------------------------------------------------------------
!  keskitahtiaika Greenwichissa 0 UTC
!-------------------------------------------------------------
double precision function sidt0(jd)
  double precision, intent(in) :: jd
  double precision :: p, j

  j=jd/36525.0
  p=24.0*dmod(0.27905727D+0 + j*(100.00213904 + &
	      j*(0.000001078 + j*7.176e-11)), 1.0D+0)
  sidt0 = p
end function

!-------------------------------------------------------------
!  keskitahtiaika
!-------------------------------------------------------------
double precision function sidt(jd, h, longitude)
  double precision, intent(in) :: jd, h, longitude
  double precision :: st

  st=sidt0(jd)
  st=dmod(st + 1.002737908D+0*(h-zone) + longitude/15.0, 24.0D+0)
  if (st < 0.0) st = st+24.0
  sidt = st
end function

!---------------------------------------------------------------------
! juliaaninen paiva, y=vuosi, m=kuukausi, d=paiva
!---------------------------------------------------------------------
double precision function julian0(y,m,d)
  implicit none
  integer, intent(in) :: y, m, d

  julian0=367*y-7*(y+(m+9)/12)/4-3*((y+(m-9)/7)/100+1)/4+275*m/9 &
         +d+1721028.5D+0

end function

!-------------------------------------------------------------
! modifioitu juliaaninen paiva
! 0 = tammikuun 1.5, 2000
!-------------------------------------------------------------
double precision function julian (y, m, d)
  implicit none
  integer, intent(in) :: y, m, d

  julian=367*y-7*(y+(m+9)/12)/4-3*((y+(m-9)/7)/100+1)/4+275*m/9 &
         +d-730516.5

end function

!---------------------------------------------------------------------
! paivan numero, y=vuosi, m=kuukausi, d=paiva
! (1.1.1900 = 1)
!---------------------------------------------------------------------
integer function day_number(y,m,d)
  implicit none
  integer, intent(in) :: y, m, d

  day_number=367*(y-1900)-7*(y+(m+9)/12)/4 &
            -3*((y+(m-9)/7)/100+1)/4+275*m/9 &
            +d+3309
  return
end function

!---------------------------------------------------------------------
! Muunnetaan paivan numero (vuoden 1900 alusta) muotoon
! y=vuosi, m=kuukausi, d=paiva ja lasketaan viikonpaiva:
! 1:maanantai, 2:tiistai, ..., 6:lauantai, 0:sunnuntai
! Muokattu algoritmista: CACM, Aug 1963, vol 6, nr. 8
!---------------------------------------------------------------------
subroutine gregorian(number,y,m,d,weekday)
  implicit none
  integer, intent(in)::  number    ! paivan numero
  integer, intent(out)::  y,m,d, & ! vuosik, kuukausi, paiva
                          weekday  ! viikonpaiva
  integer :: j

  weekday=mod(number,7)
  j=number+693901
  y=(4*j-1)/146097
  j=mod(4*j-1,146097)
  d=j/4
  j=(4*d+3)/1461
  d=mod(4*d+3,1461)/4+1
  m=(5*d-3)/153
  d=mod(5*d-3,153)/5+1
  y=100*y+j
  if (m < 10) then
     m=m+3
  else
     m=m-9
     y=y+1
  endif
  return
end subroutine

!-------------------------------------------------------------
! Auringon ja Kuun paikka
! Algoritmi perustuu artikkeliin
!   T.C. van Flandern, K.F. Pulkkinen: Low-Precision Formulae
!   for Planetary Positions,
!   Astrophysical Journal Supplement, vol 41, #2, 391-411,
!   November 1979
!-------------------------------------------------------------

!-------------------------------------------------------------
!  Aurinko                                                    
!-------------------------------------------------------------
subroutine sun(T, RA, DEC, dist)
  double precision :: T,  RA, DEC, dist

  double precision :: U,V,W 

  V=0.39785*sin(C(7))                          &
   -0.01000*sin(C(7)-C(8))                     &
   +0.00333*sin(C(7)+C(8))                     &
   -0.00021*sin(C(7))*T                        &
   +0.00004*sin(C(7)+2.0*C(8))                 &
   -0.00004*cos(C(7))                          &
   -0.00004*sin(C(5)-C(7))                     &
   +0.00003*sin(C(7)-C(8))*T   
  U=1.0000                                     &
   -0.03349*cos(C(8))                          &
   -0.00014*cos(C(8)*2.0)                      &
   +0.00008*cos(C(8))*T                        &
   -0.00003*sin(C(8)-C(19))
  W=0.03211*sin(C(8))                          &
   -0.04129*sin(C(7)*2.0)                      &
   +0.00104*sin(2.0*C(7)-C(8))                 &
   -0.00035*sin(2.0*C(7)+C(8))-0.0001          &
   -0.00008*sin(C(8))*T                        &
   -0.00008*sin(C(5))                          &
   +0.00007*sin(2.0*C(8))                      &
   +0.00005*sin(2.0*C(7))*T                    &
   +0.00003*sin(C(1)-C(7))                     &
   -0.00002*cos(C(8)-C(19))                    &
   +0.00002*sin(4.0*C(8)-8.0*C(16)+3.0*C(19))  &
   -0.00002*sin(C(8)-C(13))                    &
   -0.00002*cos(2.0*C(8)-2.0*C(13))            
 
  call coords(0,V,U,W,C(7),RA,DEC)
  dist = 1.00021 * sqrt(U)

end subroutine

!-------------------------------------------------------------
! Kuu                                                         
!-------------------------------------------------------------
subroutine moon(T, RA, DEC, dist)
  double precision :: T,  RA, DEC, dist

  double precision :: U,V,W 

  V=0.39558*sin(C(3)+C(5))                         &
   +0.08200*sin(C(3))                              &
   +0.03257*sin(C(2)-C(3)-C(5))                    &
   +0.01092*sin(C(2)+C(3)+C(5))                    &
   +0.00666*sin(C(2)-C(3))                         &
   -0.00644*sin(C(2)+C(3)-2.0*C(4)+C(5))           &
   -0.00331*sin(C(3)-2.0*C(4)+C(5))                &
   -0.00304*sin(C(3)-2.0*C(4))                     &
   -0.00240*sin(C(2)-C(3)-2.0*C(4)-C(5))           &
   +0.00226*sin(C(2)+C(3))                         &
   -0.00108*sin(C(2)+C(3)-2.0*C(4))                &
   -0.00079*sin(C(3)-C(5))                         &
   +0.00078*sin(C(3)+2.0*C(4)+C(5))                &
   +0.00066*sin(C(3)+C(5)-C(8))                    &
   -0.00062*sin(C(3)+C(5)+C(8))                    &
   -0.00050*sin(C(2)-C(3)-2.0*C(4))                &
   +0.00045*sin(2.0*C(2)+C(3)+C(5))                &
   -0.00031*sin(2.0*C(2)+C(3)-2.0*C(4)+C(5))       &
   -0.00027*sin(C(2)+C(3)-2.0*C(4)+C(5))           &
   -0.00024*sin(C(3)-2.0*C(4)+C(5)+C(8))           &
   -0.00021*sin(C(3)+C(5))*T                       &
   +0.00018*sin(C(3)-C(4)+C(5))                    &
   +0.00016*sin(C(3)+2.0*C(4))                     &
   +0.00016*sin(C(2)-C(3)-C(5)-C(8))               &
   -0.00016*sin(2.0*C(2)-C(3)-C(5))                &
   -0.00015*sin(C(3)-2.0*C(4)+C(8))                &
   -0.00012*sin(C(2)-C(3)-2.0*C(4)-C(5)+C(8))      &
   -0.00011*sin(C(2)-C(3)-C(5)+C(8)) 

  U=1.00000                                        &
   -0.10828*cos(C(2))                              &
   -0.01880*cos(C(2)-2.0*C(4))                     &
   -0.01479*cos(2.0*C(4))                          &
   +0.00181*cos(2.0*C(2)-2.0*C(4))                 &
   -0.00147*cos(2.0*C(2))                          &
   -0.00105*cos(2.0*C(4)-C(8))                     &
   -0.00075*cos(C(2)-2.0*C(4)+C(8))                &
   -0.00067*cos(C(2)-C(8))                         &
   +0.00057*cos(C(4))                              &
   +0.00055*cos(C(1)+C(8))                         &
   -0.00046*cos(C(2)+2.0*C(4))                     &
   +0.00041*cos(C(2)-2.0*C(3))                     &
   +0.00024*cos(C(8))                              &
   +0.00017*cos(2.0*C(4)+C(8))                     &
   +0.00013*cos(C(2)-2.0*C(4)-C(8))                &
   -0.00010*cos(C(2)-4.0*C(4))                     

  W=0.10478*sin(C(2))                              &
   -0.04105*sin(2.0*C(3)+2.0*C(5))                 &
   -0.02130*sin(C(2)-2.0*C(4))                     &
   -0.01779*sin(2.0*C(3)+C(5))                     &
   +0.01774*sin(C(5))                              &
   +0.00987*sin(2.0*C(4))                          &
   -0.00338*sin(C(2)-2.0*C(3)-2.0*C(5))            &
   -0.00309*sin(C(8))                              &
   -0.00190*sin(2.0*C(3))                          &
   -0.00144*sin(C(2)+C(5))                         &
   -0.00144*sin(C(2)-2.0*C(3)-C(5))                &
   -0.00113*sin(C(2)+2.0*C(3)+2.0*C(5))            &
   -0.00094*sin(C(2)-2.0*C(4)+C(8))                &
   -0.00092*sin(2.0*C(2)-2.0*C(4))                 &
   +0.00071*sin(2.0*C(4)-C(8))                     &
   +0.00070*sin(2.0*C(2))                          &
   +0.00067*sin(C(2)+2.0*C(3)-2.0*C(4)+2.0*C(5))   &
   +0.00066*sin(2.0*C(3)-2.0*C(4)+C(5))            &
   -0.00066*sin(2.0*C(4)+C(5))                     &
   +0.00061*sin(C(2)-C(8))                         &
   -0.00058*sin(C(4))                              &
   -0.00049*sin(C(2)+2.0*C(3)+C(5))                &
   -0.00049*sin(C(2)-C(5))                         &
   -0.00042*sin(C(2)+C(8))                         &
   +0.00034*sin(2.0*C(3)-2.0*C(4)+2.0*C(5))        &
   -0.00026*sin(2.0*C(3)-2.0*C(4))                 &
   +0.00025*sin(C(2)-2.0*C(3)-2.0*C(4)-2.0*C(5))   &
   +0.00024*sin(C(2)-2.0*C(3))                     &
   +0.00023*sin(C(2)+2.0*C(3)-2.0*C(4)+C(5))       &
   +0.00023*sin(C(2)-2.0*C(4)-C(5))                &
   +0.00019*sin(C(2)+2.0*C(4))                     &
   +0.00012*sin(C(2)-2.0*C(4)-C(8))                &
   +0.00011*sin(C(2)-2.0*C(4)+C(5))                &
   +0.00011*sin(C(2)-2.0*C(3)-2.0*C(4)-C(5))       &
   -0.00010*sin(2.0*C(4)+C(8))

  call coords(1,V,U,W,C(1),RA,DEC)
  dist = 60.40974 * 6378.140 * sqrt(U)
end subroutine

!-------------------------------------------------------------
! apufunktioita                                               
!-------------------------------------------------------------
double precision function A(B, C, D)
  double precision, intent(in) :: B, C, D
  double precision :: E,sig,i
  E=B+C*D
  sig=1.0
  if (E < 0.0) sig = -1.0
  A = sig * dmod(abs(E), 1.0D+0)*two_pi
end function

!-------------------------------------------------------------
! INITS alustaa taulukon C annetulle paivamaaralle            
! JS=modifioitu juliaaninen pvm, JD-2451545.0                 
!-------------------------------------------------------------
subroutine inits(JS)
  double precision, intent(in) :: JS

  C(1) =A(0.606434D+0, 0.03660110129D+0,JS) 
  C(2) =A(0.374897D+0, 0.03629164709D+0,JS) 
  C(3) =A(0.259091D+0, 0.03674819520D+0,JS) 
  C(4) =A(0.827362D+0, 0.03386319198D+0, JS) 
  C(5) =A(0.347343D+0,-0.00014709391D+0, JS)
  C(7) =A(0.779072D+0, 0.00273790931D+0, JS) 
  C(8) =A(0.993126D+0, 0.00273777850D+0, JS) 
  C(9) =A(0.700695D+0, 0.01136771400D+0, JS) 
  C(10)=A(0.485541D+0, 0.01136759566D+0, JS)
  C(11)=A(0.566441D+0, 0.01136762384D+0, JS)
  C(12)=A(0.505498D+0, 0.00445046867D+0, JS)
  C(13)=A(0.140023D+0, 0.00445036173D+0, JS)
  C(14)=A(0.292498D+0, 0.00445040017D+0, JS)
  C(15)=A(0.987573D+0, 0.00145575328D+0, JS)
  C(16)=A(0.053856D+0, 0.00145561327D+0, JS)
  C(17)=A(0.849694D+0, 0.00145569456D+0, JS)
  C(18)=A(0.089608D+0, 0.00023080893D+0, JS)
  C(19)=A(0.056531D+0, 0.00023080893D+0, JS)
  C(20)=A(0.814794D+0, 0.00023080893D+0, JS)
  C(21)=A(0.133295D+0, 0.00009294371D+0, JS)
  C(22)=A(0.882987D+0, 0.00009294371D+0, JS)
  C(23)=A(0.821218D+0, 0.00009294371D+0, JS)
  C(24)=A(0.870169D+0, 0.00003269438D+0, JS)
  C(25)=A(0.400589D+0, 0.00003269438D+0, JS)
  C(26)=A(0.664614D+0, 0.00003265562D+0, JS)
  C(27)=A(0.846912D+0, 0.00001672092D+0, JS)
  C(28)=A(0.725368D+0, 0.00001672092D+0, JS)
  C(29)=A(0.480856D+0, 0.00001663715D+0, JS)
  C(31)=A(0.663854D+0, 0.00001115482D+0, JS)
  C(32)=A(0.041020D+0, 0.00001104864D+0, JS)
  C(33)=A(0.357355D+0, 0.00001104864D+0,JS)

end subroutine

!-------------------------------------------------------------
! COORDS muuntaa U,V,W ekvatoriaalisiksi koordinateiksi       
!-------------------------------------------------------------
subroutine coords(OBJ, V, U, W, XL, RA, dec)
  integer, intent(in) :: OBJ 
  double precision :: V, U, W, XL, RA, dec
  double precision :: x,y,z

  x=V ; y=U; z=W 
  W=W/sqrt(U-V*V)
  V=V/sqrt(U)

  ! this should never happen
  if (abs(W) > 1.0) then
    write(*,'(" ===> W/sqrt(U-V*V)="F10.4)') W
    if (W<0.0) then
        W=-1.0  
    else 
        W=1.0 
    end if
  end if 
  if (abs(V) > 1.0) then
    write(*,'(" ===>V/sqrt(U)=",F10.4)') V
    if (V<0.0) then 
       V= -1.0  
    else 
       V=1.0 
    end if
  end if

  RA = XL + asin(W)
  dec = asin(V)
  if (RA < 0.0) RA = RA+two_pi
  RA = RA*24.0/two_pi
  dec = dec*360.0/two_pi

end subroutine

end module


!-------------------------------------------------------------
! moduulia kayttava testiohjelma, talleta eri tiedostoon      
!-------------------------------------------------------------
program astrotest
use astro

double precision :: jd, jj, T, st, ra, dec, dist
integer::  year=2003, month=1, day=1, wd, dn, hour=0

  ! alustetaan muuttujat; tama on toistettava erikseen
  ! jokaiselle havaintohetkelle ja -paivalle
  jd = julian (year, month, day) + hour/24.0D+0
  call INITS  (jd)
  T=jd/36525.0 + 1.0 

  ! tarkistetaan viikonpaiva
  dn = day_number (year, month, day) 
  call gregorian(dn, year, month, day, wd)
  write(*,'(F10.1,I8,I3,".",I2,".",I4,I3)') &
         jd, dn, day, month, year, wd

  ! lasketaan koordinatit
  call SUN (T, ra, dec, dist)
  write(*, &
  '("RA= ",F8.5,"h, dec= ",F8.4," deg, distance=",F8.4," AU")') & 
    ra, dec, dist

  call MOON (T, ra, dec, dist)
  write(*, &
  '("RA= ",F8.5,"h, dec= ",F8.4," deg, distance=",F8.0," km")') & 
    ra, dec, dist
 
end program