program pousuite_soleil
c       calcule les ephemerides du soleil
c       a partir de la date donnee (arret par cntrl <c>)
        real*8 jour,mois,annee,heure,minutes,nj,j,h8,n,increm,latitude
        real*8 longitude
        real*16 l0,m,l1,p0,p1,e,a,k0,xx,yy,l,ld,v,p,t,u,r,rd,l2,ts,sd,de,sr
        integer h,mm,s
        pi=3.14159265359
        unit_astro=1.496e8
        latitude=45*pi/180
        longitude=0.3
        type*,'***************************************************'
        type*,'*                                                 *'
        type*,'*                   PROGRAMME                     *'
        TYPE*,'*              POURSUITE   DU SOLEIL              *'
        type*,'*                                                 *'
        type*,'*               DE 10 MIN EN 10 MIN               *'
        type*,'***************************************************'
        TYPE*,'                                                   '
c       entrees des donnees
        type*,'ENTREZ LA DATE DU DEBUT (JJ MM AAAA) :'
        read*,jour,mois,annee
        type*,'HEURE TU (HH MM)  :'
        read*,heure,minutes
c       calcul nombre jours depuis 1/1/1900
150     h8=(heure+minutes/60)/24
        j=jour+h8
        n=annee*365+31*(mois-1)+j
        if (mois.gt.2) go to 160
        annee=annee-1
160     n=n+int(annee/4)-int(annee/100)+int(annee/400)
        if (mois.le.2) go to 190
        n=n-int((mois-1)*0.4+2.7)
190     nj=n-694325    !nj=nombre de jours depuis 1/1/1901
        type*,'nombre de jours depuis le 1/1/1901  :',nj
210     l0=4.8689
        l1=0.0172027914
        p0=4.9085
        p1=8.1856e-07
        e=0.01675104
        k=7
        a=1.00000023
        p=p0+p1*nj
        l2=l0+l1*nj
        m=l2-p
        m2=m
c       equation de kepler
        u=m
        do 200 kk=0,k
        u=m+e*sin(u)
200     continue
        v=2*atan(tan(u/2)*sqrt((1+e)/(1-e)))
        r=a*(1-e*cos(u))              !r=rayon vecteur du soleil
        type*,'rayon vecteur du soleil  :',r
        type*,'soit :',r*unit_astro,' km'
        l=v+p
        xs=r*cos(l)
        ys=r*sin(l)
        ld=l*180/pi
        ld=(ld/360-int(ld/360))*360
        if (ld.lt.0) ld=ld+360
c       ld=int(ld*10+0.5)/10          !ld=longitude du soleil
C       type*,'longitude du soleil      :',ld
c       calcul temps sideral
        rd=1.7273+1.72027914e-2*nj+h8*2*pi-longitude*pi/180
        rd=(rd/2/pi-int(rd/2/pi))*2*pi
        if(rd.lt.0) rd=rd+2*pi
        ts=rd
c       calcul  asc droite et declinaison
        b=0
        ep=0.40927971
        sd=cos(ep)*sin(b)+sin(ep)*cos(b)*sin(l)
        de=atan(sd/sqrt(1-sd*sd))    !de=declinaison
        sr=-sin(ep)*sin(b)+cos(ep)*cos(b)*sin(l)
        rd=atan(sr/cos(b)/cos(l))
        if (cos(b)*cos(l).lt.0) rd=rd+pi
        if (rd.lt.0) rd=rd+(2*pi)
        ar=rd    !ar=asc droite
        h=int(rd/pi*12)
        mm=int((rd-h*pi/12)*720/pi)
        s=int((rd-h*pi/12-m*pi/720)*43200/pi)
        type*,'ascension droite         :',h,mm,s
        type*,'declinaison              :',int(de*18000/pi+0.5)/100
c       calcul azimut elevation
        ang_hor=ts-ar
        cos_dis_zen=sin(latitude)*sin(de)+cos(latitude)*cos(de)*cos(ang_hor)
        dis_zen=acos(cos_dis_zen)
        elev=(90*pi/180)-dis_zen
        elev=elev*180/pi
        type*,'elevation :',elev
        as=cos(de)*sin(ang_hor)/cos(dis_zen)
        ac=(-cos(latitude)*sin(de)+sin(latitude)*cos(de)*cos(ang_hor))/
     1cos(dis_zen)
        azim=atan(as/ac)
        if (ac.le.0)azim=azim+pi
        azim=azim*180/pi
        azim=(azim/360-int(azim/360))*360
        if(azim.lt.0)azim=azim+360
        type*,'azimut    :',azim
c       incrementation
        minutes=minutes+10
        if(minutes.ge.60)heure=heure+1
        if(minutes.ge.60)minutes=minutes-60
        if(heure.ge.24) goto 300
        type*,'***********************************************************'
        type*,'                                                           '
        go to 150
300     type*,'FIN DE LA JOURNEE'
        end