c-------------------------------------------------------------------------- c preparation de l'integration de l'equation de la precession: c creation des fichiers binaires c - pour les elements t,k,h,q,p c - pour leurs derivees calculees a l'aide de la methode c des differences centrales. c c les parametres sont pris dans le fichier NAMELIST : 'prepa.par' c nomascpos: le nom du fichier ASCII pour les temps positifs c nomascneg: " " " " " " " " ngatifs c nomfich : le nom du fichier binaire construit avec les elements c nomficder: le nom du fichier binaire cosntruit avec les derivees c datedebut: la date initiale en millions d'annees c datefin : la date finale " " " c ( -20 <= datedebut < datefin <= +10 ) c statut : 'new' ou 'unknown' c c (c) Astronomie et Systemes Dynamiques /Bureau Des Longitudes (1993) c version 0.82 (7/7/93) c J. Laskar, F. Joutel c--------------------------------------------------------------------------- subroutine DERIVEES c****************************************************************** c c SOUS PROGRAMME DU CALCUL DES DERIVEES c A PARTIR DES VALEURS TABULEES DES FONCTIONS ak,ah,aq,ap c PAR LA METHODE DES DIFFERENCES CENTRALES c c (c) Astronomie et Systemes Dynamiques /BDL (1993) c******************************************************************* implicit double precision (a-h,o-z) parameter (ndi=20,nbf=4) parameter (nbel=5,nacd=100) dimension Taux(nbel,nacd),DTaux(nbel,nacd) dimension dy(nbel) character*50 nomfich,nomfichder,nomascpos,nomascneg character *20 statut common/CSTD/datedebut,datefin,nomfich,nomfichder, &nomascpos,nomascneg,statut c----------------------------------------------------------------- c Fichiers de travail c----------------------------------------------------------------- nrecl=8*nacd*nbel open(10,file=nomfich,status='old',access='direct', *recl=nrecl,err=98,iostat=ios10) open(13,file=nomfichder,status=statut,access='direct', *recl=nrecl,err=99,iostat=ios13) c c------------------------------------------------------------------ c Donnees et Initialisation c------------------------------------------------------------------ c c initialisation du pas: xh=1000.d0 c c intervalle de travail: tpsinit=datedebut*1.d6 tpsfinal=datefin*1.d6 nbpc=dabs((tpsfinal-tpsinit)/xh/nacd) c c constantes de la methode des differences centrales: call calcon(iord,0) c c c------------------------------------------------------------------ c Evaluation des derivees c------------------------------------------------------------------ c c Afin de posseder suffisamment de donnees pour les calculs tpsinit=tpsinit-(nacd-ndi)*1.d3 tpsfinal=tpsfinal+(nacd-ndi)*1.d3 c c Evaluation pour le premier enregistrement: tps=tpsinit do i=1,ndi do j=1,nbel DTaux(j,i)=0 enddo enddo do i=ndi+1,nacd call evderiv(tps,xh,nbf,dy) DTaux(1,i)=tps*1.d-3 do j=1,nbf DTaux(j+1,i)=dy(j) enddo tps=tps+xh enddo nrec=int((tps-tpsinit)*1.d-3/nacd) +1 write(13,rec=nrec)((Dtaux(j,i),j=1,nbel),i=1,nacd) c do iu=1,nbpc do i=1,nacd call evderiv(tps,xh,nbf,dy) DTaux(1,i)=tps*1.d-3 do j=1,nbf DTaux(j+1,i)=dy(j) enddo tps=tps+xh enddo nrec=int((tps-tpsinit)*1.d-3/nacd) +1 write(13,rec=nrec)((Dtaux(j,i),j=1,nbel),i=1,nacd) enddo c c Evaluation pour le dernier enregistrement: do i=1,nacd-ndi call evderiv(tps,xh,nbf,dy) DTaux(1,i)=tps*1.d-3 do j=1,nbf DTaux(j+1,i)=dy(j) enddo tps=tps+xh enddo do i=nacd-ndi+1,nacd do j=1,nbel DTaux(j,i)=0 enddo enddo nrec=int((tps-tpsinit)*1.d-3/nacd) +1 write(13,rec=nrec)((Dtaux(j,i),j=1,nbel),i=1,nacd) goto 100 c 99 write(*,*) 'erreur ',ios13,'dans le fichier :',nomfichder stop 98 write(*,*) 'erreur ',ios10,'dans le fichier :',nomfich stop 100 write(*,*) 'The file ',nomfichder,' is created.' close(10) close(13) return end SUBROUTINE evderiv(tps,xh,nbf,dy) c******************************************************************* c c EVALUATION DE DERIVEES PAR LA METHODE DES DIFFERENCES CENTRALES c c A L'ENTREE: c tps:date a laquelle les derivees sont evaluees(en annees) c tab:tableau des elements dont on calcule les derivees c xh:pas de temps des donnees c nbf:nombre d'elements c c EN SORTIE: c dy(nbf):tableau des derivees c (c) Astronomie et Systemes Dynamiques/BdL (1993) c******************************************************************* implicit double precision (a-h,o-z) parameter (ndi=20,nbfmax=4) dimension a(-ndi:ndi),da(-ndi:ndi,0:ndi) dimension ta(2*ndi+1),ya(nbfmax,2*ndi+1) dimension dy(nbf) dimension Tel(nbfmax) dimension der(20) common/cda/a,da c eps=1.d-16 c c------------------------------------------------------------------- c Tableaux A de la methode des differences centrales c------------------------------------------------------------------- c call Cherchel(tps,ta,ya) c c------------------------------------------------------------------- c Methode des differences centrales c------------------------------------------------------------------- do 80 i=1,nbf do 70 ii=-ndi,ndi a(ii)=ya(i,ndi+ii+1) 70 continue call diff(eps,iord,0) call deriv(iord,xh,der,0) dy(i)=der(1) 80 continue return 99 end program prepa c-------------------------------------------------------------------------- c preparation de l'integration de l'equation de la precession: c creation des fichiers binaires c - pour les elements t,k,h,q,p c - pour leurs derivees calculees a l'aide de la methode c des differences centrales. c c les parametres sont pris dans le fichier NAMELIST : 'prepa.par' c nomascpos: le nom du fichier ASCII pour les temps positifs c nomascneg: " " " " " " " " n‚gatifs c nomfich : le nom du fichier binaire construit avec les elements c nomficder: le nom du fichier binaire cosntruit avec les derivees c datedebut: la date initiale en millions d'annees c datefin : la date finale " " " c ( -20 <= datedebut < datefin <= +10 ) c statut : 'new' ou 'unknown' c c (c) Astronomie et Systemes Dynamiques /BDL (1993) c------------------------------------------------------------------------- implicit double precision (a-h,p-y) character *50 nomfich,nomfichder,nomascpos,nomascneg character *20 statut common/CSTD/datedebut,datefin,nomfich,nomfichder, &nomascpos,nomascneg,statut NAMELIST/NAMSTD/nomfich,nomfichder,nomascpos,nomascneg,datedebut, &datefin,statut open(8,file='prepa.par',form='formatted',status='old') read(8,NAMSTD) close(8) write(*,*) ' Solution La90-La93 for the Earth ' write(*,*) write(*,*) ' La93.prepa version 0.8' write(*,*) ' preparation step ' write(*,*) ' (c) ASD/BdL (1993) ' write(*,*) write(*,*) ' ASCII file for positive time :',nomascpos write(*,*) ' ASCII file for negative time :',nomascneg write(*,*) ' Binary file for elliptical elements :',nomfich write(*,*) ' Binary file for derivatives of ell. el :',nomfichder write(*,*) ' starting time (Myr) :',datedebut write(*,*) ' ending time (Myr) :',datefin call TRANSBIN call DERIVEES write(*,*) ' The preparation step is completed normally' end subroutine TRANSBIN c**************************************************************************** c c Ce sous programme lit les fichiers sequentiels ASCII et les reecrit en un c fichier a acces direct binaire, avec une longueur d'enregistrement c correspondant aux donnees (t,k,h,q,p) pour 100 000 ans. c (c) Astronomie et Systemes Dynamiques /BDL (1993) c**************************************************************************** implicit double precision (a-h,o-z) parameter(nacd=100,nbel=5) dimension Taux(nbel,nacd) character *50 nomfich,nomfichder,nomascpos,nomascneg character *20 statut common/CSTD/datedebut,datefin,nomfich,nomfichder, &nomascpos,nomascneg,statut c c Les fichiers ASCII: c-------------------- open(11,file=nomascneg,status='old') open(12,file=nomascpos,status='old') c c Le fichier a acces direct: c--------------------------- nrecl=8*nbel*nacd open(10,file=nomfich,status=statut,access='direct', *recl=nrecl,err=99,iostat=ios) c c c Il faut lire suffisamment de donnees : idebut=int(datedebut)*1000.d0-nacd ifin=int(datefin)*1000.d0+nacd c c--------------------------------------------------------------------------- c Lecture et reecriture c--------------------------------------------------------------------------- c if (idebut.lt.0.and.ifin.lt.0) then ndebut=idebut/nacd nfin=ifin/nacd-1 pdebut=0 pfin=-1 endif if (idebut.gt.0.and.ifin.gt.0) then ndebut=1 nfin=0 pdebut=idebut/nacd pfin=ifin/nacd-1 endif if (idebut.lt.0.and.ifin.gt.0) then ndebut=idebut/nacd nfin=-1 pdebut=0 pfin=ifin/nacd-1 endif c c Pour les temps negatifs: c------------------------- read(11,*) do i=-1,ifin,-1 read(11,*) enddo do n=nfin,ndebut,-1 do i=nacd,1,-1 read(11,*) (Taux(j,i),j=1,nbel) enddo nrec=(n-ndebut)+1 write(10,rec=nrec)((Taux(j,i),j=1,nbel),i=1,nacd) enddo c c Pour les temps positifs: c------------------------- do i=0,idebut-1 read(12,*) enddo do n=pdebut,pfin do i=1,nacd read(12,*) (Taux(j,i),j=1,nbel) enddo nrec=(n-ndebut)+1 if (idebut.gt.0) nrec=n-pdebut+1 write(10,rec=nrec)((Taux(j,i),j=1,nbel),i=1,nacd) enddo c c goto 100 99 write(*,*) 'error ',ios,' in the file :',nomfich stop 100 write(*,*) 'The file ', nomfich, ' is created.' close(10) close(11) close(12) return end SUBROUTINE Cherchel(tps,ta,ya) c******************************************************************** c c LES VALEURS t,k,h,q,p DESTINEES AUX DIFFERENCES c CENTRALES AUTOUR DE LA DATE TPS c c ENTREE: c tps:la date autour de laquelle on desire les elements (en annees) c c SORTIE: c ta(2*na+1) : tableau des temps equidistribues autour de tps, c allant de tps-na*1000 a tps+na*1000 ans. c ya(nbf,2*na+1): tableau des elements (k,h,q,p) correspondant c c (c) Astronomie et Systemes Dynamiques/BdL (1993) c******************************************************************** implicit double precision (a-h,o-z) parameter (na=20,nbel1=5) parameter (nbf=4,nacd=100) dimension ta(2*na+1),ya(nbf,2*na+1),y(nbf) dimension Taux(2,nbel1,nacd) logical first data first/.true./ data nreccour1/-1/ data nreccour2/-1/ save first,itdebut,ncour1,ncour2 save Taux,nreccour1,nreccour2 if (first .eqv. .true.) then read(10,rec=1) debut itdebut=int(debut) first=.false. endif c c Numero d'ordre de l'enregistrement qui contient tps c ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ tt=tps/1000.d0 itt=int(tt) nrec=(itt-itdebut)/nacd+1 itdeb=itdebut+(nrec-1)*nacd iplace=itt-itdeb+1 if (iplace.lt.na+1) then idebord=-1 else if (iplace.gt.nacd-na) then idebord=+1 else idebord=0 endif endif c c Pas de debordement c ^^^^^^^^^^^^^^^^^^^ if (idebord.eq.0) then m=2*na+1 nl1=iplace-na if (nrec.eq.nreccour1) then n1=ncour1 else if (nrec.eq.nreccour2) then n1=ncour2 else n1=1 ncour1=1 nreccour1=nrec call Listableau(n1,nrec,Taux,Dtaux) endif endif else nrecp=nrec+idebord if (idebord.lt.0) then nl1=nacd+(iplace-na) m=na-iplace+1 else nl1=iplace-na m=nacd-nl1+1 endif nl2=1 c c Chargement eventuel des tableaux en memoire c ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ if (nreccour1.eq.nrec) then if (nreccour2.ne.nrecp) then call calncour(ncour1,ncour2) nreccour2=nrecp call Listableau(ncour2,nrecp,Taux) endif else if (nreccour2.eq.nrec) then if (nreccour1.ne.nrecp) then call calncour(ncour2,ncour1) nreccour1=nrecp call Listableau(ncour1,nrecp,Taux) endif else if (nreccour1.eq.nrecp) then call calncour(ncour1,ncour2) nreccour2=nrec call Listableau(ncour2,nrec,Taux) else if (nreccour2.eq.nrecp) then call calncour(ncour2,ncour1) nreccour1=nrec call Listableau(ncour1,nrec,Taux) else ncour1=1 ncour2=2 nreccour1=nrec nreccour2=nrecp call Listableau(ncour1,nrec,Taux) call Listableau(ncour2,nrecp,Taux) endif endif endif endif c c recherche du premier tableau c ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ if (nreccour1.eq.nrec) then if (idebord.lt.0) then n1=ncour2 n2=ncour1 else n1=ncour1 n2=ncour2 endif else if (idebord.lt.0) then n1=ncour1 n2=ncour2 else n1=ncour2 n2=ncour1 endif endif endif c c chargement du tableau d'interpolation c ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 100 call Remplistableau(n1,nl1,m,n2,nl2,taux,ta,ya) return end subroutine Listableau(n,nrec,Taux) c----------------------------------------------------------------------------- c Lis la partie 'n' du tableau Taux c----------------------------------------------------------------------------- implicit double precision (a-h,p-y) parameter (nacd=100,nbel1=5) dimension Taux(2,nbel1,nacd) character *50 nomfich,nomfichder,nomascpos,nomascneg character*20 statut common/CSTD/datedebut,datefin,nomfich,nomfichder, &nomascpos,nomascneg,statut c read(10,rec=nrec,err=99,iostat=ios) * ((Taux(n,i,j),i=1,nbel1),j=1,nacd) go to 100 c 99 write(*,*) 'Erreur ',ios,'dans le fichier ',nomfich stop c 100 return end subroutine Remplistableau(n1,nl1,m,n2,nl2,Taux,ta,ya) c----------------------------------------------------------------------------- c Emplit les tableaux ta et ya a l'aide de Taux(n1, , ) sur une longueur m c a partir de l'indice nl1 puis de Taux(n2, , ) sur une longueur 2*na+1-m a c partir de l'indice nl2 c----------------------------------------------------------------------------- implicit double precision (a-h,o-z) parameter (na=20,nbel1=5,nbf=4,nacd=100) dimension ta(2*na+1),ya(nbf,2*na+1) dimension Taux(2,nbel1,nacd) dimension DTaux(2,nbel1,nacd) c nl1p=nl1-1 nl2p=nl2-1 do i=1,m ta(i)=Taux(n1,1,nl1p+i)*1.D3 do j=1,nbel1-1 ya(j,i)=Taux(n1,j+1,nl1p+i) end do end do if (m.lt.2*na+1) then do i=1,2*na+1-m ta(m+i)=Taux(n2,1,nl2p+i)*1.D3 do j=1,nbel1-1 ya(j,m+i)=Taux(n2,j+1,nl2p+i) end do end do endif return end subroutine calncour(n,m) c----------------------------------------------------------------------------- c m est l'indice du tableau libre (m=3-n) c----------------------------------------------------------------------------- implicit double precision (a-h,p-y) if (n.eq.1) then m=2 else m=1 endif return end SUBROUTINE DIFF(EPS,IORD,IPRT) ********************************************************************** * * calcule le tableau des differences centrales d'un tableau A(-NA:NA) * * A(-NA:NA) (In) Tableau des valeurs tabulees d'une fonction * EPS (In) valeurs relative a partir de laquelle les differences * n'ont plus de sens (typiquement EPS=1.D-13) * DA(-NA:NA,0:NA) (In/Out) tableau des differences centrales * IORD (Out) ordre maximum des diferences calculees obtenu a * partir de EPS * en sortie, on a donc uniquement DA(-NA:NA,0:IORD) * IPRT (In) impression si IPRT=1 * * Remarque delta^{2k+1} y_{j+1/2}=DA(j,2k+1) * delta^{2k} y_{j} =DA(j,2k) * * (c) Astronomie et Systemes Dynamiques/BdL (1993) *********************************************************************** IMPLICIT DOUBLE PRECISION (A-H,O-Y) PARAMETER(NA=20) dimension A(-NA:NA),DA(-NA:NA,0:NA) DIMENSION TN(20,0:10) COMMON/CONSTN/TN COMMON/CPRINT/NFICH COMMON/CDA/A,DA SAVE IMIN=-NA IMAX=NA KMAX=NA IORD=NA ****************! calcul de la valeur moyenne des A(i) AMOY=0.D0 DO 10 I=IMIN,IMAX AMOY=AMOY+ABS(A(I)) 10 CONTINUE AMOY=AMOY/(IMAX-IMIN+1) ****************! initialisation de DA DO 20 I=IMIN,IMAX DA(I,0)=A(I) DO 30 K=1,KMAX DA(I,K)=0.D0 30 CONTINUE 20 CONTINUE ****************! calcul de DA DO 40 K=1,KMAX EPSK=EPS*2**(K-1) ITEST=0 IF (K.GT.IORD) GOTO 999 IF (MOD(K,2).EQ.0) THEN IMIN=IMIN+1 DO 50 I=IMIN,IMAX DA(I,K)=DA(I,K-1)-DA(I-1,K-1) IF (ABS(DA(I,K))/AMOY.LE.EPSK) THEN ITEST=ITEST+1 IF(ITEST.GE.3) THEN IORD=K ENDIF ENDIF 50 CONTINUE ELSE IMAX=IMAX-1 DO 55 I=IMIN,IMAX DA(I,K)=DA(I+1,K-1)-DA(I,K-1) IF (ABS(DA(I,K))/AMOY.LE.EPSK) THEN ITEST=ITEST+1 IF(ITEST.GE.3) THEN IORD=K ENDIF ENDIF 55 CONTINUE ENDIF 40 CONTINUE 999 CONTINUE ****************! impression de DA si IPRT=1 IF (IPRT.EQ.1) THEN WRITE(NFICH,*) 'Tableau des differences centrales' WRITE(NFICH,*) ' IORD=',IORD Write(NFICH,1000) NFOIS=(IORD+1)/8+1 DO 60 N=1,NFOIS KD=(N-1)*8 KF=MIN(N*8-1,IORD) WRITE(NFICH,1002) (K,K=KD,KF) DO 70 I=-NA,NA WRITE(NFICH,1001) I,(DA(I,K),K=KD,KF) 70 CONTINUE Write(NFICH,1000) 60 CONTINUE ENDIF 1000 FORMAT (1X,///,1X) 1001 FORMAT (1X,I4,8D15.6) 1002 FORMAT (1X,4X,8(6X,I3,6X)) END SUBROUTINE CALCON(IORD,IPRT) ************************************************************************ * * calcule les constante la methode des differences centrales * (H. MINEUR p59) * * N_k^{k+2j}=TN(k,j) * ************************************************************************ IMPLICIT DOUBLE PRECISION(A-H,O-Y) dimension A(-20:20),DA(-20:20,0:20) DIMENSION TN(20,0:10) COMMON/CONSTN/TN COMMON/CPRINT/NFICH COMMON/CDA/A,DA SAVE DO 10 K=1,20 TN(K,0)=1.D0 DO 20 j=1,10 TN(K,J)=0.D0 20 CONTINUE 10 CONTINUE DO 30 J=1,10 TN(1,J)=-(2.D0*J-1)**2/(8.D0*J*(2*J+1))*TN(1,J-1) TN(2,J)=-J**2/(2.D0*(J+1)*(2*J+1))*TN(2,J-1) 30 CONTINUE DO 40 K=3,20 DO 50 J=1,10 H=K+2*J TN(K,J)=(K*(K-1)*TN(K-2,J) $ -0.25D0*(H-2)**2*TN(K,J-1)) /(H*(H-1)) 50 CONTINUE 40 CONTINUE ******************! IMPRESSION SI IPRT=1 IF (IPRT.EQ.1) THEN NFOIS=(10+1)/4+1 DO 70 N=1,NFOIS JD=4*(N-1) JF=JD+4-1 DO 60 K=1,20 WRITE(NFICH,*) (K,(K+2*J),TN(K,J),J=JD,JF) 60 CONTINUE 70 CONTINUE ENDIF 1000 FORMAT (1X,4(2I3,D25.16)) END SUBROUTINE DERIV(IORD,XH,DER,IPRT) c*********************************************************************** IMPLICIT DOUBLE PRECISION (A-H,O-Y) PARAMETER(NA=20) dimension A(-NA:NA),DA(-NA:NA,0:NA) DIMENSION TN(20,0:10),DER(20) COMMON/CONSTN/TN COMMON/CPRINT/NFICH COMMON/CDA/A,DA SAVE RPI=ATAN2(1.D0,0.D0)*2.D0 DO 10 K=1,20 DER(K)=0.D0 10 CONTINUE DO 20 K=1,20 JMAX=(IORD-K)/2 IF (JMAX.LE.0) GOTO 999 IF (IPRT.EQ.1) THEN WRITE(NFICH,*) ENDIF IF (MOD(K,2).EQ.0) THEN DO 30 J=0,JMAX DER(K)=DER(K)+TN(K,J)*DA(0,K+2*J) IF (IPRT.EQ.1) THEN WRITE(NFICH,1000) K,J,DER(K),DER(K)/XH**K ENDIF 30 CONTINUE ELSE DO 40 J=0,JMAX DER(K)=DER(K)+(K+1+2*J)/(1.D0*(K+1))*TN(K+1,J) & *(DA(0,K+2*J)+DA(-1,K+2*J))/2.D0 IF (IPRT.EQ.1) THEN WRITE(NFICH,1000) K,J,DER(K),DER(K)/XH**K ENDIF 40 CONTINUE ENDIF DER(K)=DER(K)/XH**K 999 CONTINUE 20 CONTINUE 1000 FORMAT(1X,2I5,2F30.18) END