program rayon_vect_sol
c       calcule le rayon vecteur du soleil
c       a la date dite
        real*8 jour,mois,annee,heure,minutes,nj,h8,n,j
        real*16 l0,l1,p0,p1,e,a,k0,ray,p,l2,m,u
        unit_astro=1.496e8
        type*,'CALCUL DU RAYON VECTEUR DU SOLEIL'
        type*,'ENTREZ LA DATE (HH MM AAAA) :'
        read*,jour,mois,annee
        type*,'HEURE TU  (HH MM) :'
        read*,heure,minutes
        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
        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
200     u=m+e*sin(u)
        ray=a*(1-e*cos(u))
        type*,'rayon vecteur du soleil :',ray,'u a'
        type*,'soit :',ray*unit_astro,' km'
        end