Program Code for Simulation of Non-Linear Plasma Oscillations
/******************************************************************************
* This Fortran Program simulates the non-linear oscillations of a magnetized
* collisionless plasma induced by electromagnetic waves. It is a slight
* modification of the original program used to obtain the results in my
* corresponding paper in as far as it has an additional iteration step
* (similar to a second order Runge-Kutta method) which considerably improves
* the accuracy for a given integration stepwidth.
*
* The resulting time series for the x and y components of the amplitude
* of the oscillation (referred to the initial Larmor circle) as well as the
* velocity vector is stored on file for further processing
*
* Input parameters are:
* MOD: if ='CONT': results of the present run are appended to existing file
* otherwise: results are written to new file
* DEN: plasma density (in cm-3
* DOM: frequency offset from the resonant frequency (in Hz)
* ELFS: electric field amplitude of the driving field (in cgs-units)
* MFS: field strength of the magnetic background field (in Gauss)
* X0: x-coordinate of the starting position (on arbitrary circle)
* Y0: y-coordinate of the starting position (on arbitrary circle)
* VX0: x-coordinate of the starting velocity (in cm/sec)
* VY0: y-coordinate of the starting velocity (in cm/sec)
* DDT: time stepwidth for integration (in sec)
* DDF: time stepwidth for storing result (in sec)
* DT: time stepwidth for this run (in sec)
* T0: rise/fall time for the driving field (in sec)
* TOUT: time at which driving field is switched off (with decay time T0) (in sec)
*
* Note: For a physically realistic result one needs to average the
* results for different starting positions (i.e. different values of
* X0,Y0 and VX0,VY0) on the Larmor circle. 4-8 equidistant positions on
* the circle have proved to be sufficient in this respect.
*
* I did some spot checks of the results obtained with this improved version
* (using the NFortran Light compiler (which itself is based on the Gnu Fortran
* G77 compiler)), but unlike for the original version, I have
* not been able to produce and analyze results on a larger scale.
*
* If there are any problems, please let me know and I mention the issue here.
* (for contact details see http://www.plasmaphysics.org.uk/feedback.htm ).
*
* Thomas Smid, June 2007
*********************************************************************************/
PROGRAM PLASMAOSC
IMPLICIT REAL*8 (A-Z)
INTEGER*4 N,NTF,NT,NTT,NTN,CTR,NOUT
CHARACTER*4 MOD
CHARACTER*3 STAT
CHARACTER*10 ACC
NAMELIST /INPUT/ MOD,DEN,DOM,ELFS,MFS,X0,Y0,VX0,VY0,
& DT,DDT,DDF,T0,TOUT
DATA PI/3.141592653589793D00/,EC/-4.8032D-10/,EM/9.1095D-28/,
& VL/2.9979D10/
OPEN (UNIT=22,FILE='INFILE')
READ(22,NML=INPUT)
IF (MOD .EQ. 'CONT') THEN
STAT='OLD'
ACC='APPEND '
ELSE
STAT='NEW'
ACC='SEQUENTIAL'
END IF
OPEN (UNIT=23,FILE='OUTFILE',STATUS=STAT,ACCESS=ACC,
&FORM='FORMATTED', RECL=18)
OPEN (UNIT=24,FILE='CONTFILE',STATUS=STAT,FORM='FORMATTED')
OPEN (UNIT=25,FILE='CTRFILE',STATUS=STAT,RECL=2,FORM='FORMATTED')
OPEN (UNIT=26,FILE='CONTPFILE',STATUS=STAT,ACCESS=ACC,
&FORM='FORMATTED',RECL=90)
NTN=IDNINT(DDF/DDT)
DDF=NTN*DDT
NT=IDNINT(DT/DDT)
c **** If MOD='CONT', check if previous program run has terminated normally, ****
c **** If yes, read its last coordinate values and use as new start coordinates ****
IF (MOD .EQ. 'CONT') THEN
READ(25,17) CTR
IF (CTR .EQ.0) THEN
GO TO 99
END IF
READ (24,14) X,Y,VX,VY,T
CLOSE (UNIT=24)
REWIND (UNIT=25)
END IF
CTR=0
WRITE(25,17) CTR
c **** Define some constants ****
PI2=PI*2
PIH=PI/2
EMC00=EC/EM
EMC0=EC/EM*ELFS
EMCVL0=EMC0/VL
OMP=DSQRT(4*PI*DEN*EC**2/EM)
OMP20=OMP**2
OMP2=OMP20
OML=EC*MFS/EM/VL
OM=DSQRT(OMP20+OML**2)-OML+DOM
LAM=2*PI*VL/OM
K=OM/VL
c **** Copy input data to file header in first program run ****
IF (MOD .NE. 'CONT') THEN
WRITE(23,16) DEN,DOM,ELFS,MFS,X0,Y0,VX0,VY0,
& DDT,DDF,T0
END IF
IF (MOD .EQ. 'CONT') THEN
GO TO 2
END IF
c **** Fix velocity sign for first program run ****
IF (X0*EC .LT. 0) THEN
VY0=DABS(VY0)
ELSE
VY0=-DABS(VY0)
END IF
IF (Y0*EC .LT. 0) THEN
VX0=-DABS(VX0)
ELSE
VX0=DABS(VX0)
END IF
VX=VX0
VY=VY0
2 CONTINUE
c **** Calculate speed and Larmor radius ****
V0=DSQRT(VX0**2+VY0**2)
RL=DABS(V0/OML)
c **** Calculate actual start coordinates for first program run ****
c **** and initialize deviation from Larmor circle ****
IF (MOD .NE. 'CONT') THEN
X=DABS(RL)*DCOS(DATAN2(Y0,X0))
Y=DABS(RL)*DSIN(DATAN2(Y0,X0))
T=0
END IF
R=DSQRT(X**2+Y**2)
RRL=1-RL/R
OMRRL=OMP2*RRL
c **** Write further file header data in first program run ****
IF (MOD .NE. 'CONT') THEN
WRITE(23,16) OM,LAM,OMP,OML,X,Y,RL,VX0,VY0,V0
END IF
c **** Initialize field- strength and phase,
c **** and x- and y- deviation from Larmor circle ****
EMF=0
EPHAS=OM*T-K*Y
DRX=X*RRL
DRY=Y*RRL
c **** Write further file header data in first program run ****
IF (MOD .NE. 'CONT') THEN
WRITE(23,16) T,EMF,DRX,DRY,VX,VY
END IF
c **** Initialize phase and amplitude of driving field ****
FSIN=DSIN(EPHAS)
FT=(1-DEXP(-T/T0))
IF (T .GT. TOUT) THEN
FT=FT*DEXP(-(T-TOUT)/T0)
END IF
EMF0=EMC0*FSIN*FT
EMF1=EMF0
NOUT=NTN
c **** Integration loop ****
DO 10 N=1,NT
AX=-OMRRL*X+EMF1+OML*VY
AY=-OMRRL*Y-VX*OML
DX1=VX*DDT+AX/2*DDT2
DY1=VY*DDT+AY/2*DDT2
DVX1=AX*DDT
DVY1=AY*DDT
X1=X+DX1
Y1=Y+DY1
VX1=VX+DVX1
VY1=VY+DVY1
T=T+DDT
EPHAS=OM*T-K*Y
FSIN=DSIN(EPHAS)
FT=(1-DEXP(-T/T0))
IF (T .GT. TOUT) THEN
FT=FT*DEXP(-(T-TOUT)/T0)
END IF
EMF2=EMC0*FSIN*FT
R=DSQRT(X1**2+Y1**2)
RRL=1-RL/R
OMRRL=OMP2*RRL
AX2=-OMRRL*X1+EMF2+OML*VY1
AY2=-OMRRL*Y1-VX1*OML
DX2=VX1*DDT+AX2/2*DDT2
DY2=VY1*DDT+AY2/2*DDT2
DVX2=AX2*DDT
DVY2=AY2*DDT
DX=(DX1+DX2)/2.0
DY=(DY1+DY2)/2.0
DVX=(DVX1+DVX2)/2.0
DVY=(DVY1+DVY2)/2.0
X=X+DX
Y=Y+DY
VX=VX+DVX
VY=VY+DVY
R=DSQRT(X**2+Y**2)
RRL=1-RL/R
OMRRL=OMP2*RRL
IF (N .EQ. NOUT) THEN
DRX=X*RRL
DRY=Y*RRL
EMF=EMF0/EMC00
WRITE(23,16) T,EMF,DRX,DRY,VX,VY
NOUT=N+NTN
EMF0=EMF2
END IF
EMF1=EMF2
10 CONTINUE
OPEN (UNIT=24,FILE='CONTFILE',STATUS='NEW')
WRITE (24,14) X,Y,VX,VY,T
WRITE (26,14) X,Y,VX,VY,T
CTR=1
REWIND(UNIT=25)
WRITE(25,17) CTR
14 FORMAT(5(1X,1PE17.10))
16 FORMAT((1X,1PE17.10))
17 FORMAT(1X,I1)
99 CLOSE (UNIT=22)
CLOSE (UNIT=23)
CLOSE (UNIT=24)
CLOSE (UNIT=25)
CLOSE (UNIT=26)
STOP
END
|