C******************************************************* C C DRIVER for ROWMBS4 (SEVEN BODY PROBLEM) C C******************************************************* C C******************************************************* C C THE RIGHTHANDSIDE OF ANDREWS SQUEEZING MECHANISM C HAS BEEN TAKEN FROM THE CODE HEM5 OF BRASEY/HAIRER C C******************************************************* IMPLICIT REAL*8 (A-H,O-Z) PARAMETER (NDIM=7,MDIM=6,NMDIM=13) PARAMETER ( LWORK=13*NDIM+17*MDIM+nmdim**2+NDIM**2) PARAMETER ( LIWORK=NMDIM+3) DIMENSION WORK(LWORK),IWORK(LIWORK) DIMENSION YX0(NDIM),YV0(NDIM),YX(NDIM),YV(NDIM) DIMENSION RHS(MDIM) DIMENSION am(ndim,ndim) REAL*4 dtime real*4 time,timearg(2) DIMENSION YXEX(7),YVEX(7),YXEXW(7),YVEXW(7) EXTERNAL FSTU,MSTU,GSTUP,GSTU C T0=0.0d0 T1 = 0.03D0 n=7 m=6 CALL DATSTU(N,M,T,YX0,YV0) NP1=N+1 NPM=N+M OPEN(8,FILE="w4") OPEN(9,FILE="w4g") write(8,*) ' ' write(9,*) ' ' close(8) close(9) YXEX( 1)= 15.810771195154 YXEX( 2)= -15.756371058412 YXEX( 3)= 4.0822240119641D-02 YXEX( 4)= -0.53473011634209 YXEX( 5)= 0.52440996587995 YXEX( 6)= 0.53473011634209 YXEX( 7)= 1.0480807410419 YVEX( 1)= 1139.9203022591 YVEX( 2)= -1424.3792951775 YVEX( 3)= 11.032911910584 YVEX( 4)= 19.293374104971 YVEX( 5)= 0.57356991482780 YVEX( 6)= -19.293374104971 YVEX( 7)= 0.32317914924787 do k=1,n yxexw(k)=dmax1(dabs(yxex(k)),1.0d0) yvexw(k)=dmax1(dabs(yvex(k)),1.0d0) end do C *************************************************** C C set the initial values for the performance test C loop C C *************************************************** tolmax=1.d-2 tolmin=1.d-5 tolstep=1.6d0 tol=tolmax C *************************************************** C C main loop C C *************************************************** 200 continue T=t0 te=t1 h0=(T1-T0)*tol**0.5d0 rtol=tol*1.d-2 atol=tol itol=1 do k=1,n YX(k)=yx0(k) yv(k)=yv0(k) end do iwork(1)=0 iwork(2)=0 iwork(3)=0 iwork(4)=0 work(1)=1.0d-9 time=dtime(timearg) C *************************************************** C C CALL THE INTEGRATOR C C ******************************************************* call rowmbs(n,m,t,te,h0,yx,yv, & rtol,atol,itol, & fstu,mstu,gstup,gstu, & work,lwork,iwork,liwork,idid) C**************************************** C end of integration loop C**************************************** time=dtime(timearg) C**************************************** C performance test C**************************************** if (idid.ne.0) goto 410 errx=0.0d0 errv=0.0d0 do k=1,n errx=errx+((YX(k)-yxex(k))/yxexw(k))**2 errv=errv+((YV(k)-yvex(k))/yvexw(k))**2 end do errx=errx**(0.5d0) errv=errv**0.5d0 CALL GSTU(M,N,T,YX,RHS) errg=0.0D0 do k=1,m errg=errg+RHS(k)**2 end do errg=errg**0.5d0 CALL GSTUP(M,N,N,T,YX,AM) errgv=0.0d0 do k=1,m RHS(k)=0.0D0 do j=1,n RHS(k)=RHS(k)+am(k,j)*YV(j) end do errgv=errgv+RHS(k)**2 end do errgv=errgv**0.5d0 open(8,file='w4',access='append') open(9,file='w4g',access='append') write(8,101) time,errx,errv,tol write(9,101) time,errg,errgv,tol close(8) close(9) 101 format(1X,4(E10.4,1X)) write(*,*) 'TOLERANCE:',tol, ' acc:',iwork(2), ' rej', & iwork(3) tol=tol/tolstep if (tol.ge.tolmin) goto 200 goto 420 410 write(*,*) 'exit with IDID=',idid,' at t=',t write(*,*) 'TOLERANCE:',tol tol=tol/tolstep if (tol.ge.tolmin) goto 200 420 continue end C******************************************************************** SUBROUTINE DATSTU(N,M,T,P,V) IMPLICIT REAL*8 (A-Z) INTEGER N,M,I DIMENSION P(N),V(N) C COMMON /FIXPAR/ XB,YB, XA,YA, XC,YC, RA,RR, DA,D, MOM, C0,L0, 1 SB,SS,SA,SD,SC, EA,E, TA,ZT,TB, FA,ZF, UB,U,UA, 2 M1,I1, M2,I2, M3,I3, M4,I4, M5,I5, M6,I6, M7,I7 C T=0.D0 C P(1)= -0.617138900142764496358948458001D-01 P(2)= 0.000000000000000000000000000000D+00 P(3)= 0.455279819163070380255912382449D+00 P(4)= 0.222668390165885884674473185609D+00 P(5)= 0.487364979543842550225598953530D+00 P(6)= -0.222668390165885884674473185609D+00 P(7)= 0.123054744454982119249735015568D+01 C DO 10 I=1,N V(I)=0.D0 10 CONTINUE CALL PARFIX RETURN END c********************************************************************** SUBROUTINE PARFIX IMPLICIT REAL*8 (A-Z) COMMON /FIXPAR/ XB,YB, XA,YA, XC,YC, RA,RR, DA,D, MOM, C0,L0, 1 SB,SS,SA,SD,SC, EA,E, TA,ZT,TB, FA,ZF, UB,U,UA, 2 M1,I1, M2,I2, M3,I3, M4,I4, M5,I5, M6,I6, M7,I7 c XB=-0.03635D0 YB= 0.03273D0 c XA=-0.06934D0 YA=-0.00227D0 c XC= 0.014D0 YC= 0.072D0 c RA= 0.00092D0 RR= 0.007D0 c DA= 0.0115D0 D = 0.028D0 c SB= 0.01043D0 SS= 0.035D0 SA= 0.01874D0 SD= 0.02D0 SC= 0.018D0 c EA= 0.01421D0 E = 0.02D0 c TA= 0.02308D0 ZT= 0.04D0 TB= 0.00916D0 c FA= 0.01421D0 ZF= 0.02D0 c UB= 0.00449D0 U = 0.04D0 UA= 0.01228D0 c M1= 0.04325D0 I1= 2.194D-6 c M2= 0.00365D0 I2= 4.41D-7 c M3= 0.02373D0 I3= 5.255D-6 c M4= 0.00706D0 I4= 5.667D-7 c M5= 0.07050D0 I5= 1.169D-5 c M6= 0.00706D0 I6= 5.667D-7 c M7= 0.05498D0 I7= 1.912D-5 c MOM = 0.033D0 c C0= 4530.D0 L0= 0.07785D0 c RETURN END C C SUBROUTINE MSTU(N,NDIM,T,P,M) IMPLICIT REAL*8 (A-Z) INTEGER N,NDIM,I,J DIMENSION P(N),M(NDIM,N) C COMMON /FIXPAR/ XB,YB, XA,YA, XC,YC, RA,RR, DA,D, MOM, C0,L0, 1 SB,SS,SA,SD,SC, EA,E, TA,ZT,TB, FA,ZF, UB,U,UA, 2 M1,I1, M2,I2, M3,I3, M4,I4, M5,I5, M6,I6, M7,I7 C TH=P(2) PH=P(4) OM=P(6) C COTH=COS(TH) SIPH=SIN(PH) SIOM=SIN(OM) C DO 10 I = 1,N DO 10 J = 1,I 10 M(I,J) = 0.D0 C M(1,1) = M1*RA**2 + M2*(RR**2-2*DA*RR*COTH+DA**2) + I1 + I2 M(2,1) = M2*(DA**2-DA*RR*COTH) + I2 M(2,2) = M2*DA**2 + I2 M(3,3) = M3*(SA**2+SB**2) + I3 M(4,4) = M4*(E-EA)**2 + I4 M(5,4) = M4*((E-EA)**2+ZT*(E-EA)*SIPH) + I4 M(5,5) = M4*(ZT**2+2*ZT*(E-EA)*SIPH+(E-EA)**2) + M5*(TA**2+TB**2) & + I4 + I5 M(6,6) = M6*(ZF-FA)**2 + I6 M(7,6) = M6*((ZF-FA)**2-U*(ZF-FA)*SIOM) + I6 M(7,7) = M6*((ZF-FA)**2-2*U*(ZF-FA)*SIOM+U**2) + M7*(UA**2+UB**2) & + I6 + I7 DO 30 J = 2,N DO 30 I = 1,J-1 30 M(I,J) = M(J,I) RETURN END C C SUBROUTINE FSTU(N,T,P,V,F) IMPLICIT REAL*8 (A-Z) INTEGER N,I,IERR DIMENSION P(N),V(N),F(N) C COMMON /FIXPAR/ XB,YB, XA,YA, XC,YC, RA,RR, DA,D, MOM, C0,L0, 1 SB,SS,SA,SD,SC, EA,E, TA,ZT,TB, FA,ZF, UB,U,UA, 2 M1,I1, M2,I2, M3,I3, M4,I4, M5,I5, M6,I6, M7,I7 C TH=P(2) GA=P(3) PH=P(4) OM=P(6) BEP=V(1) THP=V(2) PHP=V(4) DEP=V(5) OMP=V(6) EPP=V(7) c COGA=COS(GA) SIGA=SIN(GA) SITH=SIN(TH) COPH=COS(PH) COOM=COS(OM) c XD = SD*COGA + SC*SIGA + XB YD = SD*SIGA - SC*COGA + YB LANG = SQRT ((XD-XC)**2 + (YD-YC)**2) FORCE = - C0 * (LANG - L0)/LANG FX = FORCE * (XD-XC) FY = FORCE * (YD-YC) C F(1) = MOM - M2*DA*RR*THP*(THP+2*BEP)*SITH F(2) = M2*DA*RR*BEP**2*SITH F(3) = FX*(SC*COGA - SD*SIGA) + FY*(SD*COGA + SC*SIGA) F(4) = M4*ZT*(E-EA)*DEP**2*COPH F(5) = - M4*ZT*(E-EA)*PHP*(PHP+2*DEP)*COPH F(6) = - M6*U*(ZF-FA)*EPP**2*COOM F(7) = M6*U*(ZF-FA)*OMP*(OMP+2*EPP)*COOM RETURN END C C SUBROUTINE GSTU(M,N,T,P,G) IMPLICIT REAL*8 (A-Z) INTEGER M,N DIMENSION P(N),G(M) COMMON /FIXPAR/ XB,YB, XA,YA, XC,YC, RA,RR, DA,D, MOM, C0,L0, 1 SB,SS,SA,SD,SC, EA,E, TA,ZT,TB, FA,ZF, UB,U,UA, 2 M1,I1, M2,I2, M3,I3, M4,I4, M5,I5, M6,I6, M7,I7 C BE=P(1) TH=P(2) GA=P(3) PH=P(4) DE=P(5) OM=P(6) EP=P(7) c COBE=COS(BE) SIBE=SIN(BE) COGA=COS(GA) SIGA=SIN(GA) COBETH=COS(BE+TH) SIBETH=SIN(BE+TH) COPHDE=COS(PH+DE) SIPHDE=SIN(PH+DE) CODE=COS(DE) SIDE=SIN(DE) COOMEP=COS(OM+EP) SIOMEP=SIN(OM+EP) COEP=COS(EP) SIEP=SIN(EP) C G(1) = RR*COBE - D*COBETH - SS*SIGA - XB G(2) = RR*SIBE - D*SIBETH + SS*COGA - YB G(3) = RR*COBE - D*COBETH - E*SIPHDE - ZT*CODE - XA G(4) = RR*SIBE - D*SIBETH + E*COPHDE - ZT*SIDE - YA G(5) = RR*COBE - D*COBETH - ZF*COOMEP - U*SIEP - XA G(6) = RR*SIBE - D*SIBETH - ZF*SIOMEP + U*COEP - YA RETURN END C C SUBROUTINE GSTUP(M,N,MDIM,T,P,GP) IMPLICIT REAL*8 (A-Z) INTEGER I,J,M,N,MDIM DIMENSION P(N),GP(MDIM,N) COMMON /FIXPAR/ XB,YB, XA,YA, XC,YC, RA,RR, DA,D, MOM, C0,L0, 1 SB,SS,SA,SD,SC, EA,E, TA,ZT,TB, FA,ZF, UB,U,UA, 2 M1,I1, M2,I2, M3,I3, M4,I4, M5,I5, M6,I6, M7,I7 C BE=P(1) TH=P(2) GA=P(3) PH=P(4) DE=P(5) OM=P(6) EP=P(7) C COBE=COS(BE) SIBE=SIN(BE) COGA=COS(GA) SIGA=SIN(GA) COBETH=COS(BE+TH) SIBETH=SIN(BE+TH) COPHDE=COS(PH+DE) SIPHDE=SIN(PH+DE) CODE=COS(DE) SIDE=SIN(DE) COOMEP=COS(OM+EP) SIOMEP=SIN(OM+EP) COEP=COS(EP) SIEP=SIN(EP) C DO 10 I = 1,6 DO 10 J = 1,7 10 GP(I,J) = 0.D0 C GP(1,1) = - RR*SIBE + D*SIBETH GP(1,2) = D*SIBETH GP(1,3) = - SS*COGA GP(2,1) = RR*COBE - D*COBETH GP(2,2) = - D*COBETH GP(2,3) = - SS*SIGA GP(3,1) = - RR*SIBE + D*SIBETH GP(3,2) = D*SIBETH GP(3,4) = - E*COPHDE GP(3,5) = - E*COPHDE + ZT*SIDE GP(4,1) = RR*COBE - D*COBETH GP(4,2) = - D*COBETH GP(4,4) = - E*SIPHDE GP(4,5) = - E*SIPHDE - ZT*CODE GP(5,1) = - RR*SIBE + D*SIBETH GP(5,2) = D*SIBETH GP(5,6) = ZF*SIOMEP GP(5,7) = ZF*SIOMEP - U*COEP GP(6,1) = RR*COBE - D*COBETH GP(6,2) = - D*COBETH GP(6,6) = - ZF*COOMEP GP(6,7) = - ZF*COOMEP - U*SIEP RETURN END C