************************************************************************** * integsub.f subroutines to be used with integ.f * * (c) Astronomie et Systemes Dynamiques/Bureau des Longitudes (1993) * J. Laskar, F. Joutel * version.0.8 (6/6/1993) *************************************************************************** subroutine adams(n,fcn,x0,xfin,pas,y,valyp,solout) c--------------------------------------------------------------------------- c methode d'integration d'un systeme y'=f(x,y) c par la methode d'Adams du type predicteur/correcteur c ------> algorithme ordinaire <------- c Frederic Joutel (*m/4) c (c) ASD/BdL - 5/1991 - c - modifie le 26/8/92 pour des appels multiples - c-------------------------------------------------------------------------- c------ PARAMETRES c m l'ordre de la methode (m<=16) c nmax la dimension maximale du systeme c------ ENTREE ------------------------------------------------------------- c n dimension du syteme (n<=nmax) c fcn nom (externe) de la sous routine calculant c les seconds membres de l'equation c SUBROUTINE FCN(n,x,y,f) c REAL*8 x,y(n),f(n) c f(1)= ... c x0 la valeur initiale du temps c xfin la valeur finale (esperee cf SORTIE) c pas le pas c y(n) la position du syteme a l'instant x0 c valyp(n,0:m-1) le tableau de depart de la methode c (valyp(...,m-1) les seconds membres a l'instant x0) c Solout nom (externe) de la sous routine qui permet d'obtenir c les resultats intermediaires c SUBROUTINE SOLOUT(k,x,y,n) c REAL*8 x,y(n) c founit la solution y(n) a la k-ieme etape et a c l'instant x. c Ecrire une routine vide si ces valeurs ne sont pas c desirees c------ SORTIE -------------------------------------------------------------- c xfin la valeur comprise entre x0 et xfin, la plus c proche de xfin en comptant un nombre entier de pas c y(n) la position du systeme a l'instant xfin c valyp(n,0:m-1) le tableau des seconds membres a l'arrivee c (valyp(...,m-1) les seconds membres a l'instant xfin) c---------------------------------------------------------------------------- c implicit double precision (a-h,o-y) logical first external solout,fcn c parameter(m=12) parameter(nmax=51) c dimension pred(0:m-1,2),cor(0:m-1,2) dimension y(n),f(nmax),y1(nmax),valyp(n,0:m-1) c save pred save cor save first data first /.true./ c c---- preparations initiales posneg=dsign(1.D0,xfin-x0) h=dsign(abs(pas),posneg) if (first) then call initcoef(pred,cor,m) first=.false. endif nbre=int(abs((xfin-x0)/pas)) c---- boucle principale ----------- do k=1,nbre x=x0+k*h c-------- prediction --------- do i=1,n temp=0.D0 do j=0,m-1 temp=temp+pred(j,1)*valyp(i,j)/pred(j,2) end do y1(i)=y(i)+h*temp end do c-------- evaluation ---------- call fcn(n,x,y1,f) do i=1,n do j=0,m-2 valyp(i,j)=valyp(i,j+1) end do valyp(i,m-1)=f(i) end do c-------- correction --------- do i=1,n temp=0.D0 do j=0,m-1 temp=temp+cor(j,1)*valyp(i,j)/cor(j,2) end do y1(i)=y(i)+h*temp end do c-------- evaluation ----------- call fcn(n,x,y1,f) do i=1,n valyp(i,m-1)=f(i) end do do i=1,n y(i)=y1(i) end do call solout(k,x,y,n) end do c---- fin de la boucle principale xfin=x return end subroutine depart(n,fcn,x0,pas,m,y,valyp,solout) c--------------------------------------------------------------------------- c depart d'une methode d'integration a pas multiples d'un syteme c y'=f(x,y) c ------> utilisation de la routine Dopri8 <------- c Frederic Joutel (*m/4) c (c) ASD/BdL - 5/1991 - c (vf) c-------------------------------------------------------------------------- c------ PARAMETRE --------------------------------------------------------- c nmax dimension maximale du syteme c eps l'erreur locale sur une iteration de la routine c dopri8 (eps > udp) (calcule par le programme c eps=13.D0*udp, par compatibilite avec dopri8) c------ ENTREE ------------------------------------------------------------- c n dimension du syteme (n<=nmax) c fcn nom (externe) de la sous routine calculant c les seconds membres de l'equation c SUBROUTINE FCN(n,x,y,f) c REAL*8 x,y(n),f(n) c f(1)= ... c x0 la valeur initiale du temps c pas le pas c m l'ordre de la methode c y(n) la position du syteme a l'instant x0 c------ SORTIE -------------------------------------------------------------- c x0 la valeur du temps apres (m-1) pas c y(n) la position du systeme a l'instant x0 c valyp(n,0:m-1) le tableau des seconds membres a l'arrivee c (valyp(...,m-1) les seconds membres a l'instant x0) c---------------------------------------------------------------------------- implicit double precision (a-h,o-y) c parameter (nmax=51) c external fcn,solout dimension f(nmax),y(n),valyp(n,0:m-1) c *-------- calcul de l'epsilon de la machine (UDP) udp=1.D0 do while ((1.D0+udp) .GT. 1.D0) udp=udp/2.D0 end do udp=2.D0*udp eps=13.D0*udp c call fcn(n,x0,y,f) call solout(0,x0,y,n) do i=1,n valyp(i,0)=f(i) end do C do j=1,m-1 xfin=x0+pas c------ on prend la premiere estimation du pas de facon a ce que Dopri8 c------ calcule de lui-meme la valeur optimale au depart. pas1=pas pas0=pas call dopri8(n,fcn,x0,y,xfin,eps,pas0,pas1) call fcn(n,xfin,y,f) call solout(j,xfin,y,n) do i=1,n valyp(i,j)=f(i) end do x0=xfin end do return end SUBROUTINE PINT(xa,ya,n,x,nbf,y) c ******************************** c------------------------------------------------------------------------ c c A partir de fonctions tabulees (xa,ya),ce sous-programme construit c un polynome de degre (n-1) permettant de calculer la valeur de la c fonction, y, au point, x. c (cf algorithme de Neville, Numerical recipes) c c A L'ENTREE: c xa(n):tableau de dates c ya(nbf,n):valeurs des fonctions aux dates xa(n) c n:nombre de points necessaires a l'interpolation c x:date pour laquelle on effectue l'interpolation c nbf:nombre de fonctions a interpoler c c EN SORTIE: c y(nbf):tableau des valeurs interpolees a la date x c c (c) ASD/BDL (1992) c------------------------------------------------------------------------ c c Nombre de fonctions interpolees c implicit double precision (a-h,o-z) parameter (nmax=100,nbfmax=20) dimension xa(n),ya(nbf,n),c(nmax),d(nmax) dimension y(nbf),dy(nbfmax) c do 100 nvar=1,nbf ns=1 dif=dabs(x-xa(1)) do 10 i=1,n dift=dabs(x-xa(i)) if (dift.lt.dif) then ns=i dif=dift endif c(i)=ya(nvar,i) d(i)=ya(nvar,i) 10 continue y(nvar)=ya(nvar,ns) ns=ns-1 do 20 m=1,n-1 do 30 i=1,n-m ho=xa(i)-x hp=xa(i+m)-x w=c(i+1)-d(i) den=ho-hp if (den.eq.0) then write(*,*)' pause pint' pause endif den=w/den d(i)=hp*den c(i)=ho*den 30 continue if (2*ns.lt.n-m) then dy(nvar)=c(ns+1) else dy(nvar)=d(ns) ns=ns-1 endif y(nvar)=y(nvar)+dy(nvar) 20 continue 100 continue return end SUBROUTINE DOPRI8(N,FCN,X,Y,XEND,EPS,HMAX,H) C--------------------------------------------------------- C NUMERICAL SOLUTION OF A SYSTEM OF FIRST ORDER C ORDINARY DIFFERRENTIAL EQUATIONS Y'=F(X,Y). C THIS IS AN EMBEDDED RUNGE-KUTTA METHOD OF ORDER (7)8 C DUE TO DORMAND & PRINCE (WITH STEPSIZE CONTROL). C C.F. SECTION II.6 C C PAGE 435 HAIRER, NORSETT & WANNER C C INPUT PARAMETERS C ---------------- C N DIMENSION OF THE SYSTEM ( N.LE.51) C FCN NAME (EXTERNAL) OF SUBROUTINE COMPUTING THE C FIRST DERIVATIVE F(X,Y): C SUBROUTINE FCN(N,X,Y,F) C REAL*8 X,Y(N),F(N) C F(1)=.... ETC. C X INITIAL X-VALUE C XEND FINAL X-VALUE (XEND-X POSITIVE OR NEGATIVE) C Y(N) INITIAL VALUES FOR Y C EPS LOCAL TOLERANCE C HMAX MAXIMAL STEPSIZE C H INITIAL STEPSIZE GUESS C OUTPUT PARAMETERS C ----------------- C Y(N) SOLUTION AT XEND C C EXTERNAL SUBROUTINE (TO BE SUPPLIED BY THE USER) C ------------------- C SOLOUT THIS SUBROUTINE IS CALLED AFTER EVERY C STEP C SUBROUTINE SOLDOPRI(NR,X,Y,N) C REAL*8 X,Y(N) C FURNISHES THE SOLUTION Y AT THE NR-TH C GRID-POINT X (THE INITIAL VALUE IS CON- C SIDERED AS THE FIRST GRID-POINT). C SUPPLIED A DUMMY SUBROUTINE, IF THE SOLUTION C IS NOT DESIRED AT THE INTERMEDIATE POINTS. C-------------------------------------------------------------------- IMPLICIT REAL*8 (A-H,O-Z) REAL*8 K1(51),K2(51),K3(51),K4(51),K5(51),K6(51),K7(51), 1 Y(N),Y1(51) LOGICAL REJECT,FIRST EXTERNAL FCN,SOLOUT DATA FIRST /.TRUE./ SAVE FIRST COMMON/STAT/NFCN,NSTEP,NACCPT,NREJCT C-------COMMON STAT CAN BE USED FOR STATISTICS C------- NFCN NUMBER OF FUNCTION EVALUATIONS C------- NSTEP NUMBER OF COMPUTEF STEPS C------- NACCPT NUMBER OF ACCEPTED STEPS C------- NREJCT NUMBER OF REJECTED STEPS COMMON/COEFRK/C2,C3,C4,C5,C6,C7,C8,C9,C10,C11,C12,C13, &A21,A31,A32,A41,A43,A51,A53,A54,A61,A64,A65,A71,A74,A75,A76, &A81,A84,A85,A86,A87,A91,A94,A95,A96,A97,A98,A101,A104,A105,A106, &A107,A108,A109,A111,A114,A115,A116,A117,A118,A119,A1110,A121, &A124,A125,A126,A127,A128,A129,A1210,A1211,A131,A134,A135,A136, &A137,A138,A139,A1310,A1311,B1,B6,B7,B8,B9,B10,B11,B12,B13, &BH1,BH6,BH7,BH8,BH9,BH10,BH11,BH12 DATA NMAX/2000/,UROUND/2.23d-16/ C------- NMAX MAXIMAL NUMBER OF STEPS C------- UROUND SMALLEST NUMBER SATISFYING 1.D0+UROUND>1.0 C------- (TO BE ADAPTED BY THE USER) IF (FIRST) THEN CALL COEFST FIRST=.FALSE. ENDIF POSNEG=DSIGN(1.D0,XEND-X) C------- INITIAL PREPARATIONS HMAX=DABS(HMAX) H=DMIN1(DMAX1(1.D-10,DABS(H)),HMAX) H=DSIGN(H,POSNEG) EPS=DMAX1(EPS,13.D0*UROUND) REJECT=.FALSE. NACCPT=0 NREJCT=0 NFCN=0 NSTEP=0 C CALL SOLOUT(NACCPT+1,X,Y,N) C-------BASIC INTEGRATION STEP 1 IF (NSTEP.GT.NMAX.OR.X+.03D0*H.EQ.X) GOTO 79 IF((X-XEND)*POSNEG+UROUND.GT.0D0) RETURN IF((X+H-XEND)*POSNEG.GT.0.D0) H=XEND-X CALL FCN(N,X,Y,K1) 2 if (nstep.gt.nmax.or.x+0.03D0*h.eq.x) goto 79 NSTEP=NSTEP+1 C------- THE FIRST 9 STAGES DO 22 I=1,N 22 Y1(I)=Y(I)+H*A21*K1(I) CALL FCN(N,X+C2*H,Y1,K2) DO 23 I=1,N 23 Y1(I)=Y(I)+H*(A31*K1(I)+A32*K2(I)) CALL FCN(N,X+C3*H,Y1,K3) DO 24 I=1,N 24 Y1(I)=Y(I)+H*(A41*K1(I)+A43*K3(I)) CALL FCN(N,X+C4*H,Y1,K4) DO 25 I=1,N 25 Y1(I)=Y(I)+H*(A51*K1(I)+A53*K3(I)+A54*K4(I)) CALL FCN(N,X+C5*H,Y1,K5) DO 26 I=1,N 26 Y1(I)=Y(I)+H*(A61*K1(I)+A64*K4(I)+A65*K5(I)) CALL FCN(N,X+C6*H,Y1,K6) DO 27 I=1,N 27 Y1(I)=Y(I)+H*(A71*K1(I)+A74*K4(I)+A75*K5(I)+A76*K6(I)) CALL FCN(N,X+C7*H,Y1,K7) DO 28 I=1,N 28 Y1(I)=Y(I)+H*(A81*K1(I)+A84*K4(I)+A85*K5(I)+A86*K6(I)+A87*K7(I)) CALL FCN(N,X+C8*H,Y1,K2) DO 29 I=1,N 29 Y1(I)=Y(I)+H*(A91*K1(I)+A94*K4(I)+A95*K5(I)+A96*K6(I)+A97*K7(I) & +A98*K2(I)) CALL FCN(N,X+C9*H,Y1,K3) DO 30 I=1,N 30 Y1(I)=Y(I)+H*(A101*K1(I)+A104*K4(I)+A105*K5(I)+A106*K6(I) & +A107*K7(I)+A108*K2(I)+A109*K3(I)) C------- COMPUTE INTERMEDIATE SUMS TO SAVE MEMORY DO 61 I=1,N Y11S=A111*K1(I)+A114*K4(I)+A115*K5(I)+A116*K6(I)+A117*K7(I) & +A118*K2(I)+A119*K3(I) Y12S=A121*K1(I)+A124*K4(I)+A125*K5(I)+A126*K6(I)+A127*K7(I) & +A128*K2(I)+A129*K3(I) K4(I)=A131*K1(I)+A134*K4(I)+A135*K5(I)+A136*K6(I)+A137*K7(I) & +A138*K2(I)+A139*K3(I) K5(I)=B1*K1(I)+B6*K6(I)+B7*K7(I)+B8*K2(I)+B9*K3(I) K6(I)=BH1*K1(I)+BH6*K6(I)+BH7*K7(I)+BH8*K2(I)+BH9*K3(I) K2(I)=Y11S 61 K3(I)=Y12S C----- THE LAST 4 STAGES CALL FCN(N,X+C10*H,Y1,K7) DO 31 I=1,N 31 Y1(I)=Y(I)+H*(K2(I)+A1110*K7(I)) CALL FCN(N,X+C11*H,Y1,K2) XPH=X+H DO 32 I=1,N 32 Y1(I)=Y(I)+H*(K3(I)+A1210*K7(I)+A1211*K2(I)) CALL FCN(N,XPH,Y1,K3) DO 33 I=1,N 33 Y1(I)=Y(I)+H*(K4(I)+A1310*K7(I)+A1311*K2(I)) CALL FCN(N,XPH,Y1,K4) NFCN=NFCN+13 DO 35 I=1,N K5(I)=Y(I)+H*(K5(I)+B10*K7(I)+B11*K2(I)+B12*K3(I)+B13*K4(I)) 35 K6(I)=Y(I)+H*(K6(I)+BH10*K7(I)+BH11*K2(I)+BH12*K3(I)) C-----ERROR ESTIMATION ERR=0.D0 DO 41 I=1,N DENOM=DMAX1(1.D-6,DABS(K5(I)),DABS(Y(I)),2.D0*UROUND/EPS) 41 ERR=ERR+((K5(I)-K6(I))/DENOM)**2 ERR=DSQRT(ERR/DFLOAT(N)) C-----COMPUTATION OF HNEW C-----WE REQUIRE .333 <=HNEW/W<=6. FAC=DMAX1((1.D0/6.D0),DMIN1(3.D0,(ERR/EPS)**(1.D0/8.D0)/.9D0)) HNEW=H/FAC IF(ERR.GT.EPS) GOTO 51 C-----STEP IS ACCEPTED NACCPT=NACCPT+1 DO 44 I=1,N 44 Y(I)=K5(I) X=XPH CALL SOLDOPRI(NACCPT+1,X,Y,N) IF (DABS(HNEW).GT.HMAX) HNEW=POSNEG*HMAX IF (REJECT) HNEW=POSNEG*DMIN1(DABS(HNEW),DABS(H)) REJECT=.FALSE. H=HNEW GOTO 1 C-----STEP IS REJECT 51 REJECT=.TRUE. H=HNEW IF(NACCPT.GE.1) NREJCT=NREJCT+1 NFCN=NFCN-1 GOTO 2 C-----FAIL EXIT 79 WRITE(6,979)X 979 FORMAT(' EXIT OF DOPRI8 AT X=',D16.7) RETURN END C SUBROUTINE COEFST C------ THIS ROUTINE SETS THE COEFFICIENTS FOR THE DORMAND-PRINCE C------ METHOD OF ORDER 8 WITH ERROR ESTIMATOR OF ORDER 7 AND 13 STAGES IMPLICIT REAL*8 (A-H,O-Z) COMMON/COEFRK/C2,C3,C4,C5,C6,C7,C8,C9,C10,C11,C12,C13, & A21,A31,A32,A41,A43,A51,A53,A54,A61,A64,A65,A71,A74,A75,A76, & A81,A84,A85,A86,A87,A91,A94,A95,A96,A97,A98,A101,A104,A105,A106, & A107,A108,A109,A111,A114,A115,A116,A117,A118,A119,A1110,A121, & A124,A125,A126,A127,A128,A129,A1210,A1211,A131,A134,A135,A136, & A137,A138,A139,A1310,A1311,B1,B6,B7,B8,B9,B10,B11,B12,B13, & BH1,BH6,BH7,BH8,BH9,BH10,BH11,BH12 C2=1.D0/18.D0 C3=1.D0/12.D0 C4=1.D0/8.D0 C5=5.D0/16.D0 C6=3.D0/8.D0 C7=59.D0/400.D0 C8=93.D0/200.D0 C9=5490023248.D0/9719169821.D0 C10=13.D0/20.D0 C11=1201146811.D0/1299019798.D0 C12=1.D0 C13=1.D0 A21=C2 A31=1.D0/48.D0 A32=1.D0/16.D0 A41=1.D0/32.D0 A43=3.D0/32.D0 A51=5.D0/16.D0 A53=-75.D0/64.D0 A54=-A53 A61=3.D0/80.D0 A64=3.D0/16.D0 A65=3.D0/20.D0 A71=29443841.D0/614563906.D0 A74=77736538.D0/692538347.D0 A75=-28693883.D0/1125.D6 A76=23124283.D0/18.D8 A81=16016141.D0/946692911.D0 A84=61564180.D0/158732637.D0 A85=22789713.D0/633445777.D0 A86=545815736.D0/2771057229.D0 A87=-180193667.D0/1043307555.D0 A91=39632708.D0/573591083.D0 A94=-433636366.D0/683701615.D0 A95=-421739975.D0/2616292301.D0 A96=100302831.D0/723423059.D0 A97=790204164.D0/839813087.D0 A98=800635310.D0/3783071287.D0 A101=246121993.D0/1340847787.D0 A104=-37695042795.D0/15268766246.D0 A105=-309121744.D0/1061227803.D0 A106=-12992083.D0/490766935.D0 A107=6005943493.D0/2108947869.D0 A108=393006217.D0/1396673457.D0 A109=123872331.D0/1001029789.D0 A111=-1028468189.D0/846180014.D0 A114=8478235783.D0/508512852.D0 A115=1311729495.D0/1432422823.D0 A116=-10304129995.D0/1701304382.D0 A117=-48777925059.D0/3047939560.D0 A118=15336726248.D0/1032824649.D0 A119=-45442868181.D0/3398467696.D0 A1110=3065993473.D0/597172653.D0 A121=185892177.D0/718116043.D0 A124=-3185094517.D0/667107341.D0 A125=-477755414.D0/1098053517.D0 A126=-703635378.D0/230739211.D0 A127=5731566787.D0/1027545527.D0 A128=5232866602.D0/850066563.D0 A129=-4093664535.D0/808688257.D0 A1210=3962137247.D0/1805957418.D0 A1211=65686358.D0/487910083.D0 A131=403863854.D0/491063109.D0 A134=-5068492393.D0/434740067.D0 A135=-411421997.D0/543043805.D0 A136=652783627.D0/914296604.D0 A137=11173962825.D0/925320556.D0 A138=-13158990841.D0/6184727034.D0 A139=3936647629.D0/1978049680.D0 A1310=-160528059.D0/685178525.D0 A1311=248638103.D0/1413531060.D0 B1=14005451.D0/335480064.D0 B6=-59238493.D0/1068277825.D0 B7=181606767.D0/758867731.D0 B8=561292985.D0/797845732.D0 B9=-1041891430.D0/1371343529.D0 B10=760417239.D0/1151165299.D0 B11=118820643.D0/751138087.D0 B12=-528747749.D0/2220607170.D0 B13=1.D0/4.D0 BH1=13451932.D0/455176623.D0 BH6=-808719846.D0/976000145.D0 BH7=1757004468.D0/5645159321.D0 BH8=656045339.D0/265891186.D0 BH9=-3867574721.D0/1518517206.D0 BH10=465885868.D0/322736535.D0 BH11=53011238.D0/667516719.D0 BH12=2.D0/45.D0 RETURN END C SUBROUTINE SOLDOPRI(K,X,Y,N) RETURN END SUBROUTINE Telor(tps,Tel) c******************************************************************** c c LES ELEMENTS ORBITAUX k,h,q,p,dk,dh,dq,dp c POUR UNE DATE DONNEE tps. c c A L'ENTREE: c tps:la date pour laquelle on desire les elements (en annees) c c EN SORTIE: c Tel(nbf):tableau des elements orbitaux(k,h,q,p) pour la date tps c c (c) ASD/BDL 2/1992 c (vf) c******************************************************************** c implicit double precision (a-h,o-z) parameter (ni=8,nbf1=4,nbel1=5) parameter (nbf=8,nbel=9,nacd=100) dimension ta(ni),ya(nbf,ni),y(nbf) dimension Tel(nbf),Taux(2,nbel1,nacd) dimension DTaux(2,nbel1,nacd) data ittprev/-2000000/ data ttprev/-2.D9/ data nreccour1/-1/ data nreccour2/-1/ save ittprev,ncour1,ncour2,ttprev save Taux,DTaux,nreccour1,nreccour2 save ta,ya,y common/date/itdebut c c-------------------------------------------------------------------- c Numero d'ordre de l'enregistrement qui contient tps c-------------------------------------------------------------------- c tt=tps/1000.d0 if (tt.eq.ttprev) then go to 201 else ttprev=tt endif c c itt est la partie entiŠre de tt c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ itt=int(tt) if ((tt.lt.0.D0).and.(tt.lt.itt)) then itt=itt-1 endif if (itt.eq.ittprev) then go to 200 endif ittprev=itt nrec=(itt-itdebut)/nacd+1 itdeb=itdebut+(nrec-1)*nacd iplace=itt-itdeb+1 if (iplace.lt.ni/2) then idebord=-1 else if (iplace.gt.nacd-ni/2) then idebord=+1 else idebord=0 endif endif c c Pas de debordement c ^^^^^^^^^^^^^^^^^^^ if (idebord.eq.0) then m=ni nl1=iplace-ni/2+1 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-ni/2)+1 m=ni/2-iplace else nl1=iplace-ni/2+1 m=nacd-nl1+1 endif nl2=1 c c Chargement eventuel des tableaux en m‚moire c ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ if (nreccour1.eq.nrec) then if (nreccour2.ne.nrecp) then call calncour(ncour1,ncour2) nreccour2=nrecp call Listableau(ncour2,nrecp,Taux,Dtaux) endif else if (nreccour2.eq.nrec) then if (nreccour1.ne.nrecp) then call calncour(ncour2,ncour1) nreccour1=nrecp call Listableau(ncour1,nrecp,Taux,Dtaux) endif else if (nreccour1.eq.nrecp) then call calncour(ncour1,ncour2) nreccour2=nrec call Listableau(ncour2,nrec,Taux,Dtaux) else if (nreccour2.eq.nrecp) then call calncour(ncour2,ncour1) nreccour1=nrec call Listableau(ncour1,nrec,Taux,Dtaux) else ncour1=1 ncour2=2 nreccour1=nrec nreccour2=nrecp call Listableau(ncour1,nrec,Taux,Dtaux) call Listableau(ncour2,nrecp,Taux,Dtaux) 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,dtaux,ta,ya) c c Quand les tableaux necessaires sont en memoire: c^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 200 call pint(ta,ya,ni,tps,nbf,y) 201 do i=1,nbf Tel(i)=y(i) enddo return end subroutine Listableau(n,nrec,Taux,Dtaux) implicit double precision (a-h,p-y) parameter (nacd=100,nbel1=5) dimension DTaux(2,nbel1,nacd) dimension Taux(2,nbel1,nacd) c read(10,rec=nrec) ((Taux(n,i,j),i=1,nbel1),j=1,nacd) read(11,rec=nrec) ((DTaux(n,i,j),i=1,nbel1),j=1,nacd) return end subroutine Remplistableau(n1,nl1,m,n2,nl2,Taux,Dtaux,ta,ya) c c emplit les tableaux ta et ya … l'aide de (D)Taux(n1, , ) sur une longueur m c … partir de l'indice nl1 puis de (D)Taux(n2, , ) sur une longueur ni-m … partir c de l'indice nl2 c implicit double precision (a-h,o-z) parameter (ni=8,nbel1=5,nbf=8,nacd=100) dimension ta(ni),ya(nbf,ni) 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) ya(nbel1-1+j,i)=DTaux(n1,j+1,nl1p+i) end do end do if (m.lt.ni) then do i=1,ni-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) ya(nbel1-1+j,m+i)=DTaux(n2,j+1,nl2p+i) end do end do endif return end subroutine calncour(n,m) implicit double precision (a-h,p-y) if (n.eq.1) then m=2 else m=1 endif return end subroutine initcoef (pred,cor,m) ****************************************************************** * Initialisation des coefficients de la methode d'Adams PECE * * d'ordre 12 * ****************************************************************** implicit double precision (a-h,p-y) dimension pred(0:m-1,2),cor(0:m-1,2) c CORRECTEUR cor( 0,1)= .4671000000000000D+04 cor( 0,2)= .7884800000000000D+06 cor( 1,1)= -.6892878100000000D+08 cor( 1,2)= .9580032000000000D+09 cor( 2,1)= .3847093270000000D+09 cor( 2,2)= .9580032000000000D+09 cor( 3,1)= -.8706474100000000D+08 cor( 3,2)= .6386688000000000D+08 cor( 4,1)= .5012899030000000D+09 cor( 4,2)= .1596672000000000D+09 cor( 5,1)= -.9191049100000000D+08 cor( 5,2)= .1774080000000000D+08 cor( 6,1)= .1007253581000000D+10 cor( 6,2)= .1596672000000000D+09 cor( 7,1)= -.1022122330000000D+09 cor( 7,2)= .1774080000000000D+08 cor( 8,1)= .3646503700000000D+08 cor( 8,2)= .9123840000000000D+07 cor( 9,1)= -.9964241300000000D+08 cor( 9,2)= .4561920000000000D+08 cor(10,1)= .1374799219000000D+10 cor(10,2)= .9580032000000000D+09 cor(11,1)= .4777223000000000D+07 cor(11,2)= .1741824000000000D+08 c PREDICTEUR pred( 0,1)= -.4777223000000000D+07 pred( 0,2)= .1741824000000000D+08 pred( 1,1)= .3008230900000000D+08 pred( 1,2)= .9123840000000000D+07 pred( 2,1)= -.1741024827100000D+11 pred( 2,2)= .9580032000000000D+09 pred( 3,1)= .9236366290000000D+09 pred( 3,2)= .1520640000000000D+08 pred( 4,1)= -.6255517490000000D+09 pred( 4,2)= .4561920000000000D+07 pred( 5,1)= .3518392888300000D+11 pred( 5,2)= .1596672000000000D+09 pred( 6,1)= -.4129027322900000D+11 pred( 6,2)= .1596672000000000D+09 pred( 7,1)= .3568989256100000D+11 pred( 7,2)= .1596672000000000D+09 pred( 8,1)= -.1506437297300000D+11 pred( 8,2)= .1064448000000000D+09 pred( 9,1)= .1232664543700000D+11 pred( 9,2)= .1916006400000000D+09 pred(10,1)= -.6477936721000000D+10 pred(10,2)= .3193344000000000D+09 pred(11,1)= .4527766399000000D+10 pred(11,2)= .9580032000000000D+09 return end subroutine fcn(nd,tps,aki,Dki) c*********************************************************************** c EVALUATION DU SECOND MEMBRE DES EQUATIONS DE LA PRECESSION c c A L'ENTREE: c nd:dimension du systeme differentiel c tps:date a laquelle on evalue le second membre c aki(nd):variables utilisees l'integration des equations de precession c aki(1)=sin(eps)cos(psi) c aki(2)=sin(eps)sin(psi) c c EN SORTIE: c Dki(nd):tableau des second membres c c*********************************************************************** c implicit double precision (a-h,o-y) parameter (nbf=8) dimension aki(nd),Dki(nd),Tel(nbf) c c c----------------------------------------------------------------------- c Elements k,h,q,p,dk,dh,dq,dp pour la date tps c----------------------------------------------------------------------- call Telor(tps,Tel) ak=Tel(1) ah=Tel(2) aq=Tel(3) ap=Tel(4) dk=Tel(5) dh=Tel(6) dq=Tel(7) dp=Tel(8) c c----------------------------------------------------------------------- c Evaluation du second membre: c----------------------------------------------------------------------- call preces(tps,ak,ah,aq,ap,dk,dh,dq,dp,aki,Dki) return end subroutine solout(k,t,y,n) c--------------------------------------------------------------------- c Sous-programme de sortie des resultats c c Ce sous-programme fournit la solution y(n) a la k-ieme etape c et a l'instant t. c (vf) c--------------------------------------------------------------------- implicit double precision (a-h,o-z) dimension y(n) character*3 ecriture common/fichier/nechant,ecriture common/indice/ind c conversion du temps en milliers d'annees a=t*1.d-3 sineps=dsqrt(y(1)**2 +y(2)**2) eps=dasin(sineps) psi=atan2(y(2),y(1)) if (mod(ind,nechant).eq.0) then if (ecriture.eq.'oui') then write(98,1000) a,eps,psi else write(*,1000) a,eps,psi endif endif ind=ind+1 c1000 format(1x,d24.16,2d24.16) 1000 format(1x,f10.0,2d24.16) return end