* Damped Simple Pendulum in EXternal Field ( T. Tamaribuchi, 12/05/2000 ) PROGRAM dspexf PARAMETER(NDAT = 10000) IMPLICIT REAL*8(A-H,O-Z) * REAL*8 points(2,NDAT) REAL xm,ym LOGICAL loop COMMON /param/x0,v0,coeffC,coeffD,Aex,Wex DATA it,n/2*0/,IDl/3/ DATA loop/.TRUE./ DATA x,v,dt,Tend/1.5D0,0D0,1D-3,0D0/ OPEN(16,file='dspexf.dat') CALL initial(x,v,t,dt,dt2,nplt,nperex,Tend,ID) CALL output(n,x,v,t) DO WHILE(loop .AND. (Tend .LT. dt2 .OR. t .LT. Tend)) CALL update(x,v,t,dt,dt2) it = it + 1 IF(mod(it,nplt) .EQ. 0) THEN CALL pltout(ID,x,v,t) IF(mod(it,nperex) .EQ. 0) THEN t = it * dt IF(ID .EQ. 3 .OR. ID .EQ. 4) THEN n = n+1 * points(1,n) = x * points(2,n) = v C CALL output(n,x,v,t) CALL pltout3(ID,x,v,t) IF(n .GE. NDAT) loop = .FALSE. ENDIF ENDIF * get mouse info CALL gwmouse(ILB,IRB,xm,ym) IF(ILB .NE. 0 .OR. IRB .NE. 0) THEN CALL output(n,x,v,t) IF(ILB .NE. 0) THEN IDl = ID CALL inpltwin(ID) IF(ID .LT. 1 .OR. ID .GT. 4) THEN loop = .FALSE. ELSE IF(nperex .LE. 1 .AND. & (ID .EQ. 3 .OR. ID .EQ. 4)) & ID = IDl IF(ID .NE. IDl) THEN IF(ID .GE. 1 .AND. ID .LE. 4) THEN CALL resetwin(ID) n = 0 ENDIF ENDIF ENDIF ELSE IF(IRB .NE. 0) THEN CALL resetwin(ID) n = 0 DO WHILE(IRB .NE. 0) CALL gwmouse(ILB,IRB,xm,ym) ENDDO ENDIF ENDIF ENDIF ENDDO * IF(IDl .EQ. 3 .OR. IDl .EQ. 4) THEN * WRITE(16,'(I19/3F19.15/4F19.15)') * & n,dt,x0,v0,coeffC,coeffD,Aex,Wex * WRITE(16,'(I5,2F20.15)') (j,(points(i,j),i=1,2),j=1,n) * ENDIF CLOSE(16) CALL gwquit(IR) END C SUBROUTINE initial(x,v,t,dt,dt2,nplt,nperex,Tend,ID) IMPLICIT REAL*8(A-H,O-Z) REAL szMrk PARAMETER(PI = 3.14159265358979D0, NDIV = 500) COMMON /param/x0,v0,coeffC,coeffD,Aex,Wex COMMON /disp/t0,Tmax,Xmax,Vmax,VmaxL,szMrk COMMON /last/xl,vl,tl DATA NP/10/ t = 0 WRITE(*,'(/A)') & 'F = - sin(x) - C * v - D * abs(v)*v + A*cos(w*t)' WRITE(*,'(A,$)') & 'parameters: C, D, A, w = ' ** READ(*,*) coeffC,coeffD,Aex,Wex WRITE(*,*) coeffC,coeffD,Aex,Wex ! <=> ** WRITE(*,'(A,$)') & 'initial condition: x(0), v(0) = ' ** READ(*,*) x,v WRITE(*,*) x,v ! <=> ** WRITE(*,'(A,$)') 'dt, duration (in reduced unit) = ' ** READ(*,*) dt,Tend WRITE(*,*) dt,Tend ! <=> ** CALL inpltwin(ID) IF(ID .LT. 1 .OR. ID .GT. 4) ID = 1 * E0 = v*v/2 + (1-cos(x)) Tu = atan(1D0)*8 Tmax = Tu*NP Xmax = PI Vmax = sqrt(E0*2)*2 VmaxL = Vmax szMrk = real(Vmax/400) * plt_time = Tmax/NDIV IF(Wex .GT. 1D-5) THEN ext_per = 2*PI/Wex npltp = ext_per/plt_time nplt = (nint(ext_per/dt)+npltp-1)/npltp nperex = nplt*npltp dt = ext_per/nperex ELSE nplt = plt_time/dt nperex = 1 IF(ID .EQ. 3 .OR. ID .EQ. 4) ID = 1 ENDIF dt2 = dt/2 x0 = x v0 = v t0 = t xl = x vl = v tl = t * WRITE(*,'(A6,A15,2A24/)') & '#','time','position','velocity' * CALL gwopen(NW,0) CALL resetwin(ID) END C SUBROUTINE resetwin(ID) IMPLICIT REAL*8(A-H,O-Z) REAL szMrk COMMON /disp/t0,Tmax,Xmax,Vmax,VmaxL,szMrk IF(ID .EQ. 1) THEN CALL reset1 ELSE CALL reset2 ENDIF IF(ID .EQ. 4) CALL gwsetmrk(IR,1,szM,4,-1,5) END C SUBROUTINE output(n,x,v,t) IMPLICIT REAL*8(A-H,O-Z) WRITE(*,'(I6,F15.5,4X,1P2G24.15)') n,t,x,v END C SUBROUTINE update(x,v,t,dt,dt2) IMPLICIT REAL*8(A-H,O-Z) REAL*8 k1v,k1x,k2x,k2v,k3x,k3v,k4x,k4v PARAMETER(PI = 3.141592653589793D0, PI2 = 2D0*PI) k1v = accel(x,v,t)*dt k1x = v*dt k2v = accel(x+k1x/2,v+k1v/2,t+dt2)*dt k2x = (v+k1v/2)*dt k3v = accel(x+k2x/2,v+k2v/2,t+dt2)*dt k3x = (v+k2v/2)*dt k4v = accel(x+k3x,v+k3v,t+dt)*dt k4x = (v+k3v)*dt x = x + (k1x+2*k2x+2*k3x+k4x)/6 v = v + (k1v+2*k2v+2*k3v+k4v)/6 t = t + dt IF(x .GT. PI) x = x - PI2 IF(x .LE. -PI) x = x + PI2 END C REAL*8 FUNCTION accel(x,v,t) IMPLICIT REAL*8(A-H,O-Z) COMMON /param/x0,v0,C,D,A,w accel = - sin(x) - C * v - D * abs(v)*v + A*cos(w*t) END C SUBROUTINE inpltwin(ID) WRITE(*,'(A,$)') &'input: 1(x & v), 2(phase space), 3(Poincare sector), 4(2+3) = ' READ(*,*) ID END C SUBROUTINE pltout(ID,x,v,t) IMPLICIT REAL*8(A-H,O-Z) REAL szMrk COMMON /disp/t0,Tmax,Xmax,Vmax,VmaxL,szMrk COMMON /last/xl,vl,tl IF(ID .EQ. 2 .OR. ID .EQ. 4) THEN CALL pltout2(x,v,t) ELSE IF(ID .EQ. 1) THEN CALL pltout1(x,v,t) ENDIF xl = x vl = v tl = t END C SUBROUTINE pltout1(x,v,t) IMPLICIT REAL*8(A-H,O-Z) REAL szMrk PARAMETER(PI = 3.141592653589793D0) COMMON /disp/t0,Tmax,Xmax,Vmax,VmaxL,szMrk COMMON /last/xl,vl,tl IF(t .GT. t0 + Tmax) THEN t0 = t0 + Tmax CALL reset1 ENDIF CALL gwsetpen(IR,3,-1,-1,-1) IF(abs(x - xl) .LT. PI) + CALL gwline(IR,real(tl),real(xl/Xmax),real(t),real(x/Xmax)) CALL gwsetpen(IR,4,-1,-1,-1) CALL gwline(IR,real(tl),real(vl/Vmax),real(t),real(v/Vmax)) IF(abs(v) .GT. VmaxL) VmaxL = v END C SUBROUTINE reset1 IMPLICIT REAL*8(A-H,O-Z) REAL x, szMrk PARAMETER(PI = 3.14159265358979D0, PI2 = PI*2) COMMON /disp/t0,Tmax,Xmax,Vmax,VmaxL,szMrk COMMON /param/x0,v0,coeffC,coeffD,Aex,Wex Vmax = VmaxL CALL gwclear(IR,-1) CALL gwindow(IR,real(t0),-1.1,real(t0+Tmax),1.1) CALL gwsetpen(IR,0,-1,-1,-1) CALL gwline(IR,real(t0),0.0,real(t0+Tmax),0.0) n1 = Wex*t0/PI2+0.99 n2 = Wex*(t0+Tmax)/PI2 IF(n2 .GT. 0) THEN DO n = n1, n2 CALL gwline(IR,real(n*PI2/Wex),-1.1,real(n*PI2/Wex),1.1) ENDDO ENDIF i1 = int(t0 + 1)/10*10 i2 = int(t0+Tmax)/10*10 DO i = i1, i2, 10 x = i CALL gwline(IR,x,-0.05,x,0.05) CALL NUMBER(x+0.2, -0.06, 0.05, x, 0.0, -1) ENDDO END C SUBROUTINE pltout2(x,v,t) IMPLICIT REAL*8(A-H,O-Z) REAL szMrk PARAMETER(PI = 3.141592653589793D0) COMMON /disp/t0,Tmax,Xmax,Vmax,VmaxL,szMrk COMMON /last/xl,vl,tl CALL gwsetpen(IR,3,-1,-1,-1) IF(abs(x - xl) .LT. PI) + CALL gwline(IR,real(xl),real(vl),real(x),real(v)) IF(abs(v) .GT. VmaxL) VmaxL = v END C SUBROUTINE reset2 IMPLICIT REAL*8(A-H,O-Z) REAL Xwh,Vwh, szMrk COMMON /disp/t0,Tmax,Xmax,Vmax,VmaxL,szMrk Vmax = VmaxL Xwh = Xmax*1.1 Vwh = Vmax*1.1 CALL gwclear(IR,-1) CALL gwsetpen(IR,0,-1,-1,-1) CALL gwindow(IR,-Xwh, -Vwh, Xwh, Vwh) CALL gwline(IR, -Xwh, 0.0, Xwh, 0.0) CALL gwline(IR, 0.0, -Vwh, 0.0, Vwh) END C SUBROUTINE pltout3(ID,x,v,t) IMPLICIT REAL*8(A-H,O-Z) REAL szMrk COMMON /disp/t0,Tmax,Xmax,Vmax,VmaxL,szMrk IF(ID .EQ. 4) THEN CALL gwputmrk(IR,real(x),real(v)) ELSE CALL gwsetpxl(IR,real(x),real(v),4) ENDIF IF(abs(v) .GT. VmaxL) VmaxL = v END C BLOCKDATA IMPLICIT REAL*8(A-H,O-Z) COMMON /param/x0,v0,coeffC,coeffD,Aex,Wex DATA x0,v0,coeffC,coeffD,Aex,Wex + /1.5D0,0.0D0,0.25D0,0.0D0,0.7D0,0.6D0/ END