! This is the main part of the simulator of the particle in the quantum dot ! The quantum dot is modelled as a stadium with the walls described by ! the equation y**2=-A*x**2+B*x**4+C where A > 0, B > 0 , C > (A**2)/(4*B) ! __________________________________________________________________________ PROGRAM DOT IMPLICIT NONE INTEGER :: i_,j_,k_,l_ ! counters INTEGER :: SYMBOL,SYMBOL_1,SYMBOL_2,INITIAL_SYMBOL INTEGER :: num_N,num_N_max INTEGER :: N_attempts,attempt REAL(8) :: A_,B_,C_ ! parameters of the dot PARAMETER ( A_=10.0 , B_=1.0 , C_=26.0 ) REAL(8) :: pi PARAMETER (pi=3.14159265358979) REAL(8) :: X_,Y_ REAL(8) :: V_x,V_y REAL(8) :: theta REAL(8) :: VALUE CHARACTER :: ch_num TYPE :: POINT REAL(8) :: X_,Y_ REAL(8) :: phi REAL(8) :: V_x,V_y REAL(8) :: theta END TYPE POINT TYPE(POINT) :: INITIAL_POINT TYPE(POINT) :: THIS_POINT,THAT_POINT TYPE(POINT) :: BOUNCE_POINT OPEN(5,FILE='traject.data') DO i_=1,20 CALL RANDOM_INITIAL_POINT(INITIAL_POINT) THIS_POINT=INITIAL_POINT CALL JUMP(THIS_POINT,THAT_POINT,SYMBOL) BOUNCE_POINT=THAT_POINT IF(SYMBOL /=0 )THEN WRITE(5,'(F10.4" "F10.4" "F10.4" "F10.4" "I1)')& &INITIAL_POINT%X_ , INITIAL_POINT%Y_,& &BOUNCE_POINT%X_-INITIAL_POINT%X_ ,& &BOUNCE_POINT%Y_-INITIAL_POINT%Y_,SYMBOL_1 j_=1 DO WHILE(SYMBOL /= 0 .AND. j_ < 10) THIS_POINT=THAT_POINT BOUNCE_POINT=THIS_POINT CALL BOUNCE(THIS_POINT,THAT_POINT) THIS_POINT=THAT_POINT CALL JUMP(THIS_POINT,THAT_POINT,SYMBOL) WRITE(5,'(F10.4" "F10.4" "F10.4" "F10.4" "I1)')& &BOUNCE_POINT%X_ , BOUNCE_POINT%Y_,& &THAT_POINT%X_-BOUNCE_POINT%X_ ,& &THAT_POINT%Y_-BOUNCE_POINT%Y_,SYMBOL j_=j_+1 END DO WRITE(5,'(1X)') END IF END DO CLOSE(5,STATUS='keep') !__________________________________________________________________________ ! Now we are going to plot Poincare section for trajectories that ! have num_N points on the surface of the dot i.e. survived after ! num_N bounces. The initial point is the first bounce, ! points are recorded after reflection. ! The Poincare section in our case is the absolute value of coordinate X ! versus the X component of velocity at the points of collision ! recorded after the reflections. !__________________________________________________________________________ num_N_max=6 DO num_N=1,num_N_max WRITE(ch_num,'(I1)')num_N OPEN(6,FILE='poincare_'//ch_num//'.data', STATUS='REPLACE' ) CLOSE(6) END DO N_attempts=400000 ATTEMPTS: DO attempt=1,N_attempts 25 CONTINUE CALL RNDM_POINT_ON_THE_SURFACE(INITIAL_POINT) VALUE=EVALUATE(INITIAL_POINT,INITIAL_SYMBOL) THIS_POINT=INITIAL_POINT CALL JUMP(THIS_POINT,THAT_POINT,SYMBOL) IF(SYMBOL == 0 )THEN GO TO 25 END IF num_N=1 !------------------------------------------------------------------------- ! This part probably can be eliminated because poincare_1.data ! and poincare_2.data should be the same ! ( we are not interested in trajectories of the ! particles that bounced only once ) ! However I am not sure WRITE(ch_num,'(I1)')num_N OPEN(6,FILE='poincare_'//ch_num//'.data',POSITION='APPEND' ) WRITE(6,'(F10.4" "F10.4" "I1)')& &INITIAL_POINT%X_ , INITIAL_POINT%V_x , INITIAL_SYMBOL !-------------------------------------------------------------------------- DO WHILE((SYMBOL /= 0) .AND. (num_N < num_N_max)) THIS_POINT=THAT_POINT BOUNCE_POINT=THIS_POINT CALL BOUNCE(THIS_POINT,THAT_POINT) num_N=num_N+1 WRITE(ch_num,'(I1)')num_N OPEN(6,FILE='poincare_'//ch_num//'.data',POSITION='APPEND') WRITE(6,'(F10.4" "F10.4" "I1)')& &THAT_POINT%X_ , THAT_POINT%V_x , SYMBOL THIS_POINT=THAT_POINT CALL JUMP(THIS_POINT,THAT_POINT,SYMBOL) END DO END DO ATTEMPTS CONTAINS !_____________________________________________________________________! SUBROUTINE RANDOM_INITIAL_POINT(INITIAL_POINT) IMPLICIT NONE REAL(8) :: uR1_rndm,uR2_rndm,uR3_rndm REAL(8) :: sR1_rndm,sR2_rndm,sR3_rndm REAL(8) :: X_rndm,Y_rndm REAL(8) :: theta_rndm TYPE(POINT), INTENT(OUT) :: INITIAL_POINT 10 CONTINUE CALL RANDOM_SEED CALL RANDOM_NUMBER(uR1_rndm) CALL RANDOM_NUMBER(uR2_rndm) CALL RANDOM_NUMBER(uR3_rndm) sR1_rndm=2.0*(uR1_rndm-0.5) ! sR#_rndm are random real numbers sR2_rndm=2.0*(uR2_rndm-0.5) ! -1< sR#_rndm < 1 sR3_rndm=2.0*(uR3_rndm-0.5) X_rndm=4*sqrt(A_/B_)*sR1_rndm Y_rndm=4*sqrt(A_/B_)*sR2_rndm theta_rndm=2*pi*uR3_rndm IF( Y_rndm**2 < -A_*(X_rndm**2)+B_*(X_rndm**4)+C_)THEN INITIAL_POINT%X_=X_rndm INITIAL_POINT%Y_=Y_rndm INITIAL_POINT%phi=DATAN2(X_rndm,Y_rndm) INITIAL_POINT%theta=theta_rndm INITIAL_POINT%V_x=cos(theta_rndm) INITIAL_POINT%V_y=sin(theta_rndm) ELSE GO TO 10 END IF RETURN END SUBROUTINE RANDOM_INITIAL_POINT !_____________________________________________________________________! SUBROUTINE RNDM_POINT_ON_THE_SURFACE(INITIAL_POINT) ! This subroutine generates the initial point ! on the sufrace inside the dot REAL(8) :: uR1_rndm,uR2_rndm,uR3_rndm REAL(8) :: sR1_rndm,sR2_rndm,sR3_rndm REAL(8) :: X_rndm,Y_rndm REAL(8) :: phi_rndm REAL(8) :: V_x_rndm,V_y_rndm REAL(8) :: theta_rndm REAL(8) :: dFdx,t_x,t_y,n_x,n_y TYPE(POINT), INTENT(OUT) :: INITIAL_POINT CALL RANDOM_SEED CALL RANDOM_NUMBER(uR1_rndm) CALL RANDOM_NUMBER(uR2_rndm) CALL RANDOM_NUMBER(uR3_rndm) sR1_rndm=2.0*(uR1_rndm-0.5) ! sR#_rndm are random real numbers sR2_rndm=2.0*(uR2_rndm-0.5) ! -1< sR#_rndm < 1 sR3_rndm=2.0*(uR3_rndm-0.5) ! The X coordinate should be somewhere inside ! the interval -sqrt(A/2B) < X < sqrt(A/2B) ! and Y can be positive or negative, therefore ! we can write: X_rndm=SQRT(A_/(2*B_))*sR1_rndm IF(sR2_rndm > 0 )THEN Y_rndm=SQRT(-A_*(X_rndm**2)+B_*(X_rndm**4)+C_) ELSE Y_rndm=-SQRT(-A_*(X_rndm**2)+B_*(X_rndm**4)+C_) END IF ! Moreover, the velocity of the random initial point ! should be directed away from the surface , ! i.e. only "emitted" particles should be considered ! when we record components of velocities for Poincare ! section we should record them AFTER the last bounce 20 CONTINUE CALL RANDOM_NUMBER(uR3_rndm) theta_rndm=2*pi*uR3_rndm V_x_rndm=cos(theta_rndm) V_y_rndm=sin(theta_rndm) IF(Y_rndm > 0.0 )THEN dFdx=(-A_*X_rndm+2*B_*X_rndm**3)/sqrt(-A_*X_rndm**2+B_*X_rndm**4+C_) t_x=1/sqrt(1+dFdx**2) t_y=dFdx/sqrt(1+dFdx**2) n_x=t_y n_y=-t_x ELSE dFdx=-(-A_*X_rndm+2*B_*X_rndm**3)/sqrt(-A_*X_rndm**2+B_*X_rndm**4+C_) t_x=1/sqrt(1+dFdx**2) t_y=dFdx/sqrt(1+dFdx**2) n_x=-t_y n_y=t_x END IF IF(n_x*V_x_rndm+n_y*V_y_rndm .LE. 0.0 )GO TO 20 INITIAL_POINT%X_=X_rndm INITIAL_POINT%Y_=Y_rndm INITIAL_POINT%phi=DATAN2(X_rndm,Y_rndm) INITIAL_POINT%theta=theta_rndm INITIAL_POINT%V_x=V_x_rndm INITIAL_POINT%V_y=V_y_rndm RETURN END SUBROUTINE RNDM_POINT_ON_THE_SURFACE !_____________________________________________________________________! SUBROUTINE DRAW_DOT_WALLS(A_,B_,C_) IMPLICIT NONE INTEGER :: i_ INTEGER :: i_max PARAMETER (i_max=200) REAL(8), INTENT(IN) :: A_,B_,C_ REAL(8) :: X_i,dX dX=4*sqrt(A_/B_)/i_max OPEN(2,FILE='dot_walls.data') X_i=-2*sqrt(A_/B_) DO i_=1,i_max WRITE(2,'(F10.4" "F10.4)')X_i,sqrt(-A_*X_i**2+B_*X_i**4+C_) X_i=X_i+dX END DO WRITE(2,'(1X)') X_i=2*sqrt(A_/B_) DO i_=1,i_max WRITE(2,'(F10.4" "F10.4)')X_i,-sqrt(-A_*X_i**2+B_*X_i**4+C_) X_i=X_i-dX END DO CLOSE(2) RETURN END SUBROUTINE DRAW_DOT_WALLS !_____________________________________________________________________! SUBROUTINE JUMP(THE_POINT,THE_NEXT_POINT,SYMBOL) INTEGER, INTENT(OUT) :: SYMBOL REAL(8) :: delta_L1,delta_L2 REAL(8) :: delta_t,delta_tt REAL(8) :: X_,Y_ REAL(8) :: phi REAL(8) :: V_x,V_y REAL(8) :: theta REAL(8) :: R_ REAL(8) :: VALUE TYPE(POINT), INTENT(IN) :: THE_POINT TYPE(POINT), INTENT(OUT) :: THE_NEXT_POINT TYPE(POINT) :: OLD_POINT,NEW_POINT,MIN_POINT OLD_POINT=THE_POINT 425 CONTINUE CALL MOVE_POINT(OLD_POINT,NEW_POINT) IF( (NEW_POINT%X_)**2+(NEW_POINT%Y_)**2 > 16*A_/B_)& &THEN SYMBOL=0 THE_NEXT_POINT=NEW_POINT RETURN END IF 480 CONTINUE IF(ABS(NEW_POINT%V_x) > ABS(NEW_POINT%V_y))THEN CALL MIN_A(sqrt(0.02*A_/B_),NEW_POINT,MIN_POINT) CALL NEWTON_A(MIN_POINT,NEW_POINT) ELSE CALL MIN_B(sqrt(0.02*A_/B_),NEW_POINT,MIN_POINT) CALL NEWTON_B(MIN_POINT,NEW_POINT) END IF THE_NEXT_POINT=NEW_POINT VALUE=EVALUATE(THE_NEXT_POINT,SYMBOL) !print*,' abs value >',ABS(VALUE) !print*,'_______________________________________' IF( VALUE > 0.01*A_/B_ )THEN print*,'not accurate' END IF RETURN END SUBROUTINE JUMP !_____________________________________________________________________! FUNCTION F_FUNC(xi,THE_POINT) REAL(8) :: F_FUNC REAL(8), INTENT(IN) :: xi REAL(8) :: yi REAL(8) :: X_,Y_ REAL(8) :: phi REAL(8) :: V_x,V_y REAL(8) :: theta REAL(8) :: R_ REAL(8) slope TYPE(POINT), INTENT(IN) :: THE_POINT R_=sqrt(C_) X_=THE_POINT%X_ Y_=THE_POINT%Y_ phi=THE_POINT%phi V_x=THE_POINT%V_x V_y=THE_POINT%V_y theta=THE_POINT%theta IF( abs(V_x) > 0.0000001 )THEN slope=V_y/V_x ! print*,'slope=',slope ELSE WRITE(*,*)' ERROR in F_FUNC! ' WRITE(*,*)'bad slope of the trajectory' print*,'V_x=',V_x print*,'slope=',slope STOP END IF yi=Y_+slope*(xi-X_) F_FUNC=-A_*(xi**2)+B_*(xi**4)+C_-yi**2 RETURN END FUNCTION F_FUNC !_____________________________________________________________________! FUNCTION DF_FUNC(xi,THE_POINT) REAL(8) :: DF_FUNC REAL(8), INTENT(IN) :: xi REAL(8) :: yi REAL(8) :: X_,Y_ REAL(8) :: phi REAL(8) :: V_x,V_y REAL(8) :: theta REAL(8) :: R_ REAL(8) slope TYPE(POINT), INTENT(IN) :: THE_POINT R_=sqrt(C_) X_=THE_POINT%X_ Y_=THE_POINT%Y_ phi=THE_POINT%phi V_x=THE_POINT%V_x V_y=THE_POINT%V_y theta=THE_POINT%theta IF( abs(V_x) > 0.0000001 )THEN slope=V_y/V_x ! print*,'slope=',slope ELSE WRITE(*,*)' ERROR in DF_FUNC! ' WRITE(*,*)'bad slope of the trajectory' print*,'V_x=',V_x print*,'slope=',slope STOP END IF yi=Y_+slope*(xi-X_) DF_FUNC=-2*A_*xi+4*B_*(xi**3)-2*yi*slope RETURN END FUNCTION DF_FUNC !_____________________________________________________________________! FUNCTION G_FUNC(yi,THE_POINT) REAL(8) :: G_FUNC REAL(8), INTENT(IN) :: yi REAL(8) :: xi REAL(8) :: X_,Y_ REAL(8) :: phi REAL(8) :: V_x,V_y REAL(8) :: theta REAL(8) :: R_ REAL(8) slope TYPE(POINT), INTENT(IN) :: THE_POINT R_=sqrt(C_) X_=THE_POINT%X_ Y_=THE_POINT%Y_ phi=THE_POINT%phi V_x=THE_POINT%V_x V_y=THE_POINT%V_y theta=THE_POINT%theta IF(abs(V_y) > 0.0000001)THEN slope=V_x/V_y ! print*,'slope=',slope ELSE WRITE(*,*)' ERROR in G_FUNC! ' WRITE(*,*)'bad slope of the trajectory' STOP END IF xi=X_+slope*(yi-Y_) G_FUNC=-A_*(xi**2)+B_*(xi**4)+C_-yi**2 RETURN END FUNCTION G_FUNC !_____________________________________________________________________! FUNCTION DG_FUNC(yi,THE_POINT) REAL(8) :: DG_FUNC REAL(8), INTENT(IN) :: yi REAL(8) :: xi REAL(8) :: X_,Y_ REAL(8) :: phi REAL(8) :: V_x,V_y REAL(8) :: theta REAL(8) :: R_ REAL(8) slope TYPE(POINT), INTENT(IN) :: THE_POINT R_=sqrt(C_) X_=THE_POINT%X_ Y_=THE_POINT%Y_ phi=THE_POINT%phi V_x=THE_POINT%V_x V_y=THE_POINT%V_y theta=THE_POINT%theta IF(abs(V_y) > 0.0000001)THEN slope=V_x/V_y ! print*,'slope=',slope ELSE WRITE(*,*)' ERROR in DG_FUNC! ' WRITE(*,*)'bad slope of the trajectory' STOP END IF xi=X_+slope*(yi-Y_) DG_FUNC=-2*A_*xi*slope+4*B_*slope*(xi**3)-2*yi RETURN END FUNCTION DG_FUNC !_____________________________________________________________________! SUBROUTINE MIN_A(d,OLD_POINT,NEW_POINT) ! this subroutine minimizes functions of one variable ! minimization is performed with simplex method ! FUNC_F - is a funtion that is to be minimized ! x_c - is the coordinare of the central vortex ! d - is the distance between vortices of the simplex IMPLICIT NONE LOGICAL :: shift INTEGER :: N_steps_max,step PARAMETER (N_steps_max=400) REAL(8) :: x_c,d REAL(8) :: x_c_initial REAL(8) :: x_1,x_2,dd REAL(8) :: f_c,f_1,f_2 REAL(8) :: x_min,f_min TYPE(POINT), INTENT(IN) :: OLD_POINT TYPE(POINT), INTENT(OUT) :: NEW_POINT step=1 dd=d NEW_POINT=OLD_POINT x_c=OLD_POINT%X_ x_c_initial=x_c x_min=x_c 88 CONTINUE IF(ABS(F_FUNC(x_min,OLD_POINT)) > & & ABS(F_FUNC(x_c_initial,OLD_POINT)))THEN print*,'error in MIN_A!' pause END IF IF((step > N_steps_max).OR. (ABS(F_FUNC(x_min,OLD_POINT)) < 0.00001*A_/B_))THEN x_c=x_min NEW_POINT%X_=x_c NEW_POINT%Y_=OLD_POINT%Y_+(OLD_POINT%V_y/OLD_POINT%V_x)*(x_c-OLD_POINT%X_) NEW_POINT%phi=DATAN2(NEW_POINT%Y_,NEW_POINT%X_) NEW_POINT%V_x=OLD_POINT%V_x NEW_POINT%V_y=OLD_POINT%V_y NEW_POINT%theta=OLD_POINT%theta ! print 15,step,dd/d,ABS(F_FUNC(x_min,OLD_POINT)), 0.00001*A_/B_ ! 15 FORMAT(' A step=',I3,' dd/d=',G12.6,' |F|=',G12.6,' < ',G12.6) RETURN END IF shift=.FALSE. x_1=x_min-dd x_2=x_min+dd f_1=ABS(F_FUNC(x_1,OLD_POINT)) f_2=ABS(F_FUNC(x_2,OLD_POINT)) f_c=ABS(F_FUNC(x_min,OLD_POINT)) IF(f_1.LT.f_c)THEN x_min=x_1 shift=.TRUE. END IF IF(f_2.LT.f_c)THEN x_min=x_2 shift=.TRUE. END IF IF(shift)THEN x_c=x_min step=step+1 GO TO 88 ELSE dd=0.5*dd step=step+1 GO TO 88 END IF print*,'FALL THROUGH' RETURN END SUBROUTINE MIN_A !_____________________________________________________________________! SUBROUTINE MIN_B(d,OLD_POINT,NEW_POINT) ! this subroutine minimizes functions of one variable ! minimization is performed with simplex method ! ABS(FUNC_G) - is a funtion that is to be minimized ! y_c - is the coordinare of the central vortex ! d - is the distance between vortices of the simplex IMPLICIT NONE LOGICAL :: shift INTEGER :: N_steps_max,step PARAMETER (N_steps_max=400) REAL(8) :: y_c,d,y_c_initial REAL(8) :: y_1,y_2,dd REAL(8) :: f_c,f_1,f_2 REAL(8) :: y_min,f_min TYPE(POINT), INTENT(IN) :: OLD_POINT TYPE(POINT), INTENT(OUT) :: NEW_POINT NEW_POINT=OLD_POINT step=1 dd=d y_c=OLD_POINT%Y_ y_c_initial=y_c y_min=y_c 88 CONTINUE IF(ABS(G_FUNC(y_min,OLD_POINT)) >& &ABS(G_FUNC(y_c_initial,OLD_POINT)))THEN print*,'error in MIN_B!' pause END IF IF((step > N_steps_max).OR. ABS(G_FUNC(y_min,OLD_POINT)) < 0.00001*A_/B_ )THEN y_c=y_min NEW_POINT%Y_=y_c NEW_POINT%X_=OLD_POINT%X_+(OLD_POINT%V_x/OLD_POINT%V_y)*(y_c-OLD_POINT%Y_) NEW_POINT%phi=DATAN2(NEW_POINT%Y_,NEW_POINT%X_) NEW_POINT%V_x=OLD_POINT%V_x NEW_POINT%V_y=OLD_POINT%V_y NEW_POINT%theta=OLD_POINT%theta ! print 15, step, dd/d, ABS(G_FUNC(y_min,OLD_POINT)), 0.00001*A_/B_ ! 15 FORMAT(' B step=',I3,' dd/d=',G12.6,' |G|=',G12.6,' < ',G12.6 ) RETURN END IF shift=.FALSE. y_1=y_min-dd y_2=y_min+dd f_1=ABS(G_FUNC(y_1,OLD_POINT)) f_2=ABS(G_FUNC(y_2,OLD_POINT)) f_c=ABS(G_FUNC(y_min,OLD_POINT)) IF(f_1.LT.f_c)THEN y_min=y_1 shift=.TRUE. END IF IF(f_2.LT.f_c)THEN y_min=y_2 shift=.TRUE. END IF IF(shift)THEN y_c=y_min step=step+1 GO TO 88 ELSE dd=0.5*dd step=step+1 GO TO 88 END IF print*,'FALL THROUGH' RETURN END SUBROUTINE MIN_B !_____________________________________________________________________! SUBROUTINE NEWTON_A(OLD_POINT,NEW_POINT) IMPLICIT NONE REAL(8) :: X_0 TYPE(POINT), INTENT(IN) :: OLD_POINT TYPE(POINT), INTENT(OUT) :: NEW_POINT X_0=OLD_POINT%X_ IF(ABS(DF_FUNC(X_0,OLD_POINT)) < 0.000001)THEN NEW_POINT=OLD_POINT RETURN END IF NEW_POINT%X_=OLD_POINT%X_-F_FUNC(X_0,OLD_POINT)/DF_FUNC(X_0,OLD_POINT) NEW_POINT%Y_=OLD_POINT%Y_+& &(OLD_POINT%V_y/OLD_POINT%V_x)*& &(NEW_POINT%X_-OLD_POINT%X_) NEW_POINT%phi=DATAN2(NEW_POINT%Y_,NEW_POINT%X_) NEW_POINT%V_x=OLD_POINT%V_x NEW_POINT%V_y=OLD_POINT%V_y NEW_POINT%theta=OLD_POINT%theta RETURN END SUBROUTINE NEWTON_A !_____________________________________________________________________! SUBROUTINE NEWTON_B(OLD_POINT,NEW_POINT) IMPLICIT NONE REAL(8) :: Y_0 TYPE(POINT), INTENT(IN) :: OLD_POINT TYPE(POINT), INTENT(OUT) :: NEW_POINT Y_0=OLD_POINT%Y_ IF(ABS(DG_FUNC(Y_0,OLD_POINT)) < 0.000001)THEN NEW_POINT=OLD_POINT RETURN END IF NEW_POINT%Y_=OLD_POINT%Y_-G_FUNC(Y_0,OLD_POINT)/DG_FUNC(Y_0,OLD_POINT) NEW_POINT%X_=OLD_POINT%X_+& &(OLD_POINT%V_x/OLD_POINT%V_y)*& &(NEW_POINT%Y_-OLD_POINT%Y_) NEW_POINT%phi=DATAN2(NEW_POINT%Y_,NEW_POINT%X_) NEW_POINT%V_x=OLD_POINT%V_x NEW_POINT%V_y=OLD_POINT%V_y NEW_POINT%theta=OLD_POINT%theta RETURN END SUBROUTINE NEWTON_B !_____________________________________________________________________! SUBROUTINE BOUNCE(OLD_POINT,NEW_POINT) ! This subroutine gives new velocities V_x, V_y ! after bounce from the wall IMPLICIT NONE REAL(8) :: X_,Y_ REAL(8) :: V_x_old,V_x_new REAL(8) :: V_y_old,V_y_new REAL(8) :: V_n_old,V_n_new REAL(8) :: V_t_old,V_t_new REAL(8) :: theta_old,theta_new,phi REAL(8) :: n_x,n_y,t_x,t_y REAL(8) :: dFdx,DF,DX REAL(8) :: V_ TYPE(POINT), INTENT(IN) :: OLD_POINT TYPE(POINT), INTENT(OUT) :: NEW_POINT X_=OLD_POINT%X_ Y_=OLD_POINT%Y_ V_x_old=OLD_POINT%V_x V_y_old=OLD_POINT%V_y phi=OLD_POINT%phi !print*,'bounce',X_,Y_ IF(ABS(-A_*X_**2+B_*X_**4+C_-Y_**2) > 0.01*A_/B_ )THEN print*,'ERROR! Can not perform reflection.' print*,-A_*X_**2+B_*X_**4+C_-Y_**2 pause END IF IF(Y_ > 0.0 )THEN dFdx=(-A_*X_+2*B_*X_**3)/sqrt(-A_*X_**2+B_*X_**4+C_) t_x=1/sqrt(1+dFdx**2) t_y=dFdx/sqrt(1+dFdx**2) n_x=t_y n_y=-t_x ELSE dFdx=-(-A_*X_+2*B_*X_**3)/sqrt(-A_*X_**2+B_*X_**4+C_) t_x=1/sqrt(1+dFdx**2) t_y=dFdx/sqrt(1+dFdx**2) n_x=-t_y n_y=t_x END IF !OPEN(3,FILE='bounces.data',POSITION='APPEND') !WRITE(3,*)X_,Y_,n_x,n_y !CLOSE(3,status='keep') ! Now we should use the formula for mirror reflection ! Vb=Va-2(Va.n)n V_x_new=V_x_old-2*n_x*(n_x*V_x_old+n_y*V_y_old) V_y_new=V_y_old-2*n_y*(n_x*V_x_old+n_y*V_y_old) ! now let us adjust the vectors V_x_new, V_y_new ! so that the energy ( and , therefore, the absolute value of V) ! be constant V_=sqrt(V_x_new**2+V_y_new**2) V_x_new=V_x_new/V_ V_y_new=V_y_new/V_ theta_new=DATAN2(V_y_new,V_y_new) NEW_POINT%X_=X_ NEW_POINT%Y_=Y_ NEW_POINT%phi=phi NEW_POINT%V_x=V_x_new NEW_POINT%V_y=V_y_new NEW_POINT%theta=theta_new RETURN END SUBROUTINE BOUNCE FUNCTION EVALUATE(THE_POINT,SYMBOL) IMPLICIT NONE INTEGER , INTENT(OUT) :: SYMBOL REAL(8) :: EVALUATE TYPE(POINT), INTENT(IN):: THE_POINT IF(THE_POINT%X_ .GE. 0.0 .AND. THE_POINT%Y_ .GE. 0.0)SYMBOL=1 IF(THE_POINT%X_ .LT. 0.0 .AND. THE_POINT%Y_ .GE. 0.0)SYMBOL=2 IF(THE_POINT%X_ .LT. 0.0 .AND. THE_POINT%Y_ .LT. 0.0)SYMBOL=3 IF(THE_POINT%X_ .GE. 0.0 .AND. THE_POINT%Y_ .LT. 0.0)SYMBOL=4 EVALUATE=-A_*(THE_POINT%X_**2)+B_*(THE_POINT%X_**4)+C_-THE_POINT%Y_**2 RETURN END FUNCTION EVALUATE SUBROUTINE MOVE_POINT(OLD_POINT,NEW_POINT) ! This subroutine moves a given initial point OLD_POINT ! to the point NEW_POINT which is very close to the ! actual point of collision IMPLICIT NONE INTEGER :: i_step REAL(8) :: X_,Y_,X_0,Y_0,X_1,Y_1 REAL(8) :: phi REAL(8) :: V_x,V_y REAL(8) :: theta REAL(8) :: blind_step,sighted_step TYPE(POINT), INTENT(IN) ::OLD_POINT TYPE(POINT), INTENT(OUT) :: NEW_POINT X_0=OLD_POINT%X_ Y_0=OLD_POINT%Y_ phi=OLD_POINT%phi V_x=OLD_POINT%V_x V_y=OLD_POINT%V_y theta=OLD_POINT%theta blind_step=0.01*sqrt(A_/B_) sighted_step=0.005*sqrt(A_/B_) i_step=0 X_=OLD_POINT%X_ Y_=OLD_POINT%Y_ X_1=OLD_POINT%X_+V_x*(sighted_step+blind_step) Y_1=OLD_POINT%Y_+V_y*(sighted_step+blind_step) DO WHILE(-A_*(X_1**2)+B_*(X_1**4)+C_-Y_1**2 > 0.0 & &.AND.& & X_**2+Y_**2 < 16*A_/B_ & &) X_=X_0+V_x*sighted_step*i_step+V_x*blind_step Y_=Y_0+V_y*sighted_step*i_step+V_y*blind_step X_1=X_0+V_x*sighted_step*(i_step+1)+V_x*blind_step Y_1=Y_0+V_y*sighted_step*(i_step+1)+V_y*blind_step i_step=i_step+1 END DO 240 CONTINUE phi=DATAN2(X_,Y_) NEW_POINT%X_=X_ NEW_POINT%Y_=Y_ NEW_POINT%phi=phi NEW_POINT%V_x=V_x NEW_POINT%V_y=V_y NEW_POINT%theta=theta RETURN END SUBROUTINE MOVE_POINT END PROGRAM DOT