C *********************************************************
C This program calculates the refraction angle theta as well
C the shift at the window location due to the refraction as
C a function of the radial location
C *********************************************************
Parameter(npt=85,npt1=43)
Implicit Real*8 (A-H,O-Z)
Character*4 ch
character*20 datfil1,datfil2,datfil3,datfil4
Real*8 Y(4),W(4,9),rmaj,rmin,c(24),Lambda,n2,
& disprel,zk0
Real*8 rlow,deltar,rstart,zstart,deltaz,deltathe
Real*8 rarr(npt),zarr(npt),Rhoo(npt,npt),
& dRhor(npt,npt),dRhoz(npt,npt)
Real*8 R,Z,Rho,dRhodr,dRhodz,ellip,delta
Real*8 gampol1,gampol2,gampol3,tottheta,rimp
Real*4 X1,Y1,Z1
COMMON /COMBL/ rarr,zarr,Rhoo,dRhor,dRhoz
common /combl1/ alpha,omega,vellight,rmaj,rmin,
& peakdensity,zk0,pi
COMMON /COMBL2/ ellip,delta
common /prexp/ prexp
C ********************SST Parameters******************
pi = 4.0d0 * datan(1.0d0)
rmaj = 1.10d0
rmin = 0.20d0
peakdensity = 3.0d19
ellip = 1.7d0
delta = 0.6d0
C ********************SST Parameters******************
vellight = 3.0D8
Alpha = 3.17730D3 !e**2/2*Epsilon*m
C ********************Reading the profiles***************
call prof
C
c To reproduce euidensity contours below the median plane.
C
do 101 i = 1,npt
do 103 j =1,npt1-1
k1= npt-j+1
zarr(j) = -zarr(k1)
rhoo(i,j) = Rhoo(i,k1)
drhor(i,j)= drhor(i,k1)
drhoz(i,j)= -drhoz(i,k1)
103 continue
101 continue
C ********************Reading the profiles***************
freq = 6.94d11
lambda= vellight/freq
Omega = 2.0d0 * pi * freq
zk0 = omega / vellight !Free space wave vector
Print*,'lambda=', lambda
NEqn = 4
NW = NEqn
Tol = 1.0D-10
C *********************************************************
rlow = rmaj - rmin
Deltar = 1.0d-2 !increment in R
Deltaz = 1.0d-2 !increment in Z
deltathe = 0.01745d0 !increment in theta
C *********************************************************
C Prexp gives the exponential value of density profile,
C N= N0(1-RHO**2/RMIN**2)^PREXP
C *********************************************************
c Type*,'Give the exponent value for the profiles'
c Read*,prexp
nexp = 0
123 nexp= nexp+1
if (nexp .eq. 1) prexp = 0.5d0
if (nexp .eq. 2) prexp = 1.0d0
if (nexp .eq. 3) prexp = 2.0d0
if (nexp .gt. 3) stop
If(prexp .eq. 0.5d0)then
ch = 'p1b2'
elseif(prexp .eq. 1.0d0)then
ch = 'p1'
elseif(prexp .eq. 2.0d0)then
ch = 'p2'
else
ch = 'arb'
endif
c Print*,'Write <1> if it is vertical probing'
c Print*,' <2> if it is horizontal probing'
c Print*,' <3> if it is oblique probing'
c Read*,ncond
ncond = 0
1234 ncond = ncond + 1
if (ncond .gt. 3) goto 123
If(ncond .eq. 1)then
datfil1 = 'vrefd'//ch//'.dat'
open(30,File=datfil1,Status='Unknown',access='append')
elseif(ncond .eq. 2)then
datfil2 = 'hrefd'//ch//'.dat'
Open(40,File=datfil2,Status='Unknown')
elseif(ncond .eq. 3)then
datfil3 = 'orefd'//ch//'.dat'
Open(50,File=datfil3,Status='Unknown')
endif
k=0
111 Index = 1
k=k+1
If(ncond .eq. 1)then
R = rlow + dfloat(k-1)*deltar
z = -ellip*rmin
Zkrin = 0.0d0
zkzin = zk0
if (r .gt. rmin+rmaj)then
close(30)
goto 1234
endif
elseif(ncond .eq. 2 .or. ncond .eq. 3)then
k1 = 0
if (ncond .eq. 2)then
Z = dfloat(k-1) * deltaz
Theta0=0.0d0
if (z .ge. ellip*rmin)then
close(40)
goto 1234
Endif
elseif(ncond .eq. 3)then
Theta0 = dfloat(k-1)*deltathe
z = 8.85d-1*dtan(theta0)
c 8.85d-1 is distance of entry pt to outer plasma edge
if(k.gt.47)then
close(50)
goto 1234
endif
endif
r = rmaj + rmin
Zst = Z + (r-0.80d0)*dtan(theta0)
c 0.80D0 is distance of reflectors from the major axis
zkrin =- zk0*dcos(theta0)
zkzin = zk0*dsin(theta0)
endif
Y(1) = r
Y(2) = z
Y(3) = zkrin
Y(4) = zkzin
call Locate(y(1),y(2),rho,drhodr,drhodz)
c ***** Normalize k-vectors to free space vector:******
y(3) = y(3) / zk0
y(4) = y(4) / zk0
c ***** Normalize k-vectors to free space vector:******
kount = 0
Time = 0.0D0
delt = 2.0d0 * rmin/3.0D8/4000 !time step in ray trajectory
10 tend = Time + delt
kount = kount + 1
hh = delt * 0.5d0
c Use 4th-order Runge-Kutta integrator:
call rk4sys(4,time,y,delt)
c Write(*,8) (Y(i),i=1,2),zk0*y(3),zk0*y(4)
index = 2
C To locate the coordinates along the ray path
call Locate(y(1),y(2),rho,drhodr,drhodz)
If (Rho/rmin .ge. 0.995d0)then
if(ncond .eq. 1 .and. y(2) .lt. 0.0d0)ncond1=1
if(ncond .eq. 1 .and. y(2) .gt. 0.0d0)ncond1=2
if(ncond .eq. 2 .or. ncond .eq. 3)then
if(y(1) .gt. Rmaj)ncond1=1
if(y(1) .lt. Rmaj)ncond1=2
endif
If(ncond1 .eq. 1)then
rstart = y(1)
Zstart = y(2)
If(zstart .gt. rmin*ellip)then
If(ncond .eq. 2)then
close(40)
Elseif(ncond .eq. 3)then
close(50)
Endif
goto 1234
endif
Elseif(ncond1 .eq. 2)then
Y(3) = y(3) * zk0
Y(4) = y(4) * zk0
zkf = dsqrt(y(3)**2 + y(4)**2)
theta=(Zkrin*y(3)+Zkzin*y(4))/zk0/zkf
theta = dacos(theta)
Y1 = Theta
If(ncond.eq.1)then
shiftedge=Y(1)-R
Y(1)=(1.5d0-Y(2))*Y(3)/Y(4)
C Shift from the st. line trajectory without plasma
Shift=shiftedge+Y(1)
X1= R-Rmaj
Z1 = shift
If(shift .lt. 0.0d0)y1=-y1
Print*,'Vertical probing'
Print*,x1,y1,z1
Write(30,*)x1,y1,z1
Elseif(ncond.eq.2.or.ncond.eq.3)then
Shiftedge = Y(2)-zst
Y(2)=(0.80d0-Y(1))*Y(4)/Y(3)
Shift=Shiftedge+y(2)
Z1 = Shift
If(shift .lt. 0.0d0)y1=-y1
If (ncond .eq. 2)then
X1 = z
Print*,'Horizontal probing'
Print*,X1,Y1,Z1
Write(40,*)X1,Y1,Z1
Elseif (ncond .eq. 3)then
X1 = theta0
Print*,'Oblique probing'
Print*,X1*1.8d2/pi,Y1,Z1
Write(50,*)X1,Y1,Z1
Endif
Endif
goto 111
Endif
Endif
Call Gradient(NEqn,Y,dndR,dndz,Density)
ter1 = zk0 * zk0 * (Y(3)**2 + y(4)**2)*vellight**2
ter2 = Omega**2
ter3 = alpha*Density
disprel = dabs(ter1 - ter2 + ter3) / dmax1(dabs(ter1),
& dabs(ter2),dabs(ter3))
If (disprel .gt. 1.0D-3) then
Write(*,*)'Dispersion Relation is not satisfied'
endif
c if (kount.gt.10) then
c kount = 0
c Write(*,8) (Y(i),i=1,2),zk0*y(3),zk0*y(4),Density
c8 Format(' ',5(e11.4,1x))
c endif
Goto 10
End
C *******************************************************
Subroutine prof
C *******************************************************
C This subroutine calculates Rho and derivatives of Rho
C w.r.t R and z on the generated grid above the median
C plane.
C *******************************************************
Parameter(npt=85,npt1=43)
C Radial(npt) and vertical(npt1) grid points for D-shaped
C plasma above the mdeian plane
Implicit Real*8 (A-H,O-Z)
Real*8 ellip,Delta,Rmaj,Rmin,Rho,Rho1,ractual,Theta,Rx,rcal
& rdash,slope,zslope,zactual
Real*8 Rarr(npt),Zarr(npt),Rhoo(npt,npt),dRhor(npt,npt),
& dRhoz(npt,npt),zorg,rorg
logical quadrant
COMMON /COMBL/ rarr,zarr,Rhoo,dRhor,dRhoz
common /combl1/ alpha,omega,vellight,rmaj,rmin,
& peakdensity,zk0,pi
COMMON /COMBL2/ ellip,delta
common /combl3/ ractual,zactual,quadrant,zslope
External F
pi = 4.0d0 * datan(1.0d0)
Slope = - Ellip/dsin(delta)
C DeltaR and Deltaz gives the incremental grid area.
Deltar = 2.0d0*rmin/(npt-5)
Deltaz = ellip*rmin/(npt1-3)
rlow = rmaj - rmin -2*deltar
Do 100 k = 1,npt
Ractual = rlow + dfloat(k-1)*deltar
rarr(k) = ractual
Zslope = slope * ( ractual-Rmaj)
Do 101 j = 1,npt1
j2 = npt1+j-1
Zactual = dfloat(j-1) * deltaz
Zarr(j2) = zactual
Do 102 j1 = 1,3
If (j1 .eq. 2)then
Ractual = rlow + dfloat(k-1)*deltar +2.0d-4
Zactual = dfloat(j-1) * deltaz
elseif(j1 .eq. 3)then
Ractual = rlow + dfloat(k-1)*deltar
Zactual = dfloat(j-1) * deltaz + 2.0d-4
endif
103 Rholow = zactual / ellip
if (rholow.lt.1.0d-4) rholow = 1.0d-4
Rhoup = 2.0d0* rmin
Eps = 0.0d0
Nsig = 8
Maxfn = 5000
quadrant = .false.
If((zslope-zactual).gt.0.0d0.and.(ractual-rmaj)
& .lt.0.0d0)quadrant = .true.
Call zbrent(F,Eps,Nsig,Rholow,Rhoup,Maxfn,Ier)
rho = rhoup
If (dsqrt(zactual**2+(ractual-rmaj)**2).lt.1.0d-4)then
theta = pi/2
rho = zactual/ellip
else
If(zactual .gt. ellip*rho)then
Print*,'zactual=',zactual
Print*,'ellip=',ellip
Print*,'rho=',rho
endif
Theta = dasin(Zactual/ellip/Rho)
endif
If (quadrant) theta=pi-theta
Temp = theta + delta*dsin(theta)
Rcomp = rmaj + rho * dcos(temp)
zcomp = Ellip*rho*dsin(theta)
If(dabs(zcomp-zactual).gt.1.0d-4 .or.
& dabs(ractual-rcomp).gt.1.0d-4)then
c Write(30,5)ractual,rcomp,zactual,zcomp,Rho,theta*1.8d2/pi
c Write(*,*)ractual,rcomp,zactual,zcomp,Rho,theta*1.8d2/pi
c Write(30,*)'problem'
zactual= zactual+1.0d-5
goto 103
ENDIF
If (j1.eq. 1) then
Rorg = Ractual
Zorg = zactual
Rhoo(k,j2) = rho
elseif(j1 .eq. 2) then
dRhor(k,j2) = (rho - rhoo(k,j2))/2.0d-4
elseif(j1 .eq. 3)then
dRhoz(k,j2) = (rho - rhoo(k,j2))/2.0d-4
endif
if (ier.gt.128) write(*,*) ier
102 continue
c Write(40,*) rarr(k),zarr(j2),Rhoo(k,j2),drhor(k,j2),
c & drhoz(k,j2)
101 continue
100 continue
5 Format(1x,6(1x,d10.4))
end
C *************************************************************
SUBROUTINE ZBRENT (F,EPS,NSIG,A,B,MAXFN,IER)
C IMSL ROUTINE NAME - ZBRENT
C
C *************************************************************
C
C COMPUTER - IBM/SINGLE
C
C LATEST REVISION - JANUARY 1, 1978
C
C PURPOSE - ZERO OF A FUNCTION WHICH CHANGES SIGN IN A
C GIVEN INTERVAL (BRENT ALGORITHM)
C
C USAGE - CALL ZBRENT (F,EPS,NSIG,A,B,MAXFN,IER)
C
C ARGUMENTS F - AN EXTERNAL FUNCTION SUBPROGRAM F(X)
C PROVIDED BY THE USER WHICH COMPUTES F FOR
C ANY X IN THE INTERVAL (A,B). (INPUT)
C F MUST APPEAR IN AN EXTERNAL STATEMENT IN
C THE CALLING PROGRAM
C EPS - FIRST CONVERGENCE CRITERION (INPUT). A ROOT,
C B, IS ACCEPTED IF ABS(F(B)) IS LESS THAN OR
C EQUAL TO EPS. EPS MAY BE SET TO ZERO.
C NSIG - SECOND CONVERGENCE CRITERION (INPUT). A ROOT,
C B, IS ACCEPTED IF THE CURRENT APPROXIMATION
C AGREES WITH THE TRUE SOLUTION TO NSIG
C SIGNIFICANT DIGITS.
C A,B - ON INPUT, THE USER MUST SUPPLY TWO POINTS, A
C AND B, SUCH THAT F(A) AND F(B) ARE OPPOSITE
C IN SIGN.
C ON OUTPUT, BOTH A AND B ARE ALTERED. B
C WILL CONTAIN THE BEST APPROXIMATION TO THE
C ROOT OF F. SEE REMARK 1.
C MAXFN - ON INPUT, MAXFN SHOULD CONTAIN AN UPPER BOUND
C ON THE NUMBER OF FUNCTION EVALUATIONS
C REQUIRED FOR CONVERGENCE. ON OUTPUT, MAXFN
C WILL CONTAIN THE ACTUAL NUMBER OF FUNCTION
C EVALUATIONS USED.
C IER - ERROR PARAMETER. (OUTPUT)
C TERMINAL ERROR
C IER = 129 INDICATES THE ALGORITHM FAILED TO
C CONVERGE IN MAXFN EVALUATIONS.
C IER = 130 INDICATES F(A) AND F(B) HAVE THE
C SAME SIGN.
C
C PRECISION/HARDWARE - SINGLE AND DOUBLE/H32
C - SINGLE/H36,H48,H60
C
C REQD. IMSL ROUTINES - UERTST,UGETIO
C
C NOTATION - INFORMATION ON SPECIAL NOTATION AND
C CONVENTIONS IS AVAILABLE IN THE MANUAL
C INTRODUCTION OR THROUGH IMSL ROUTINE UHELP
C
C REMARKS 1. ON EXIT FROM ZBRENT, WHEN IER=0, A AND B SATISFY THE
C FOLLOWING,
C F(A)*F(B) .LE.0,
C ABS(F(B)) .LE. ABS(F(A)), AND
C EITHER ABS(F(B)) .LE. EPS OR
C ABS(A-B) .LE. MAX(ABS(B),0.1)*10.0**(-NSIG).
C THE PRESENCE OF 0.1 IN THIS ERROR CRITERION CAUSES
C LEADING ZEROES TO THE RIGHT OF THE DECIMAL POINT TO BE
C COUNTED AS SIGNIFICANT DIGITS. SCALING MAY BE REQUIRED
C IN ORDER TO ACCURATELY DETERMINE A ZERO OF SMALL
C MAGNITUDE.
C 2. ZBRENT IS GUARANTEED TO REACH CONVERGENCE WITHIN
C K = (ALOG((B-A)/D)+1.0)**2 FUNCTION EVALUATIONS WHERE
C D=MIN(OVER X IN (A,B) OF
C MAX(ABS(X),0.1)*10.0**(-NSIG)).
C THIS IS AN UPPER BOUND ON THE NUMBER OF EVALUATIONS.
C RARELY DOES THE ACTUAL NUMBER OF EVALUATIONS USED BY
C ZBRENT EXCEED SQRT(K). D CAN BE COMPUTED AS FOLLOWS,
C P = AMIN1(ABS(A),ABS(B))
C P = AMAX1(0.1,P)
C IF ((A-0.1)*(B-0.1).LT.0.0) P = 0.1
C D = P*10.0**(-NSIG)
C
C COPYRIGHT - 1977 BY IMSL, INC. ALL RIGHTS RESERVED.
C
C WARRANTY - IMSL WARRANTS ONLY THAT IMSL TESTING HAS BEEN
C APPLIED TO THIS CODE. NO OTHER WARRANTY,
C EXPRESSED OR IMPLIED, IS APPLICABLE.
C
C-----------------------------------------------------------------------
C
C SPECIFICATIONS FOR ARGUMENTS
implicit real*8 (a-h,o-z)
INTEGER NSIG,MAXFN,IER
REAL*8 F,EPS,A,B
C SPECIFICATIONS FOR LOCAL VARIABLES
INTEGER IC
REAL*8 ZERO,HALF,ONE,THREE,TEN,
1 T,FA,FB,C,FC,D,E,TOL,RM,S,P,Q,R,RONE,TEMP
DATA ZERO/0.0/,HALF/.5/,ONE/1.0/,THREE/3.0/,
1 TEN/10.0/
C FIRST EXECUTABLE STATEMENT
IER = 0
T = TEN**(-NSIG)
IC = 2
S = A
FA = F(S)
S = B
FB = F(S)
C TEST FOR SAME SIGN
IF (FA*FB.GT.ZERO) GO TO 50
5 C = A
FC = FA
D = B-C
E = D
10 IF (ABS(FC).GE.ABS(FB)) GO TO 15
A = B
B = C
C = A
FA = FB
FB = FC
FC = FA
15 CONTINUE
TOL = T*dMAX1(ABS(B),0.1d0)
RM = (C-B)*HALF
C TEST FOR FIRST CONVERGENCE CRITERIA
IF (ABS(FB).LE.EPS) GO TO 40
C TEST FOR SECOND CONVERGENCE CRITERIA
IF (ABS(C-B).LE.TOL) GO TO 40
C CHECK EVALUATION COUNTER
IF (IC.GE.MAXFN) GO TO 45
C IS BISECTION FORCED
IF (ABS(E).LT.TOL) GO TO 30
IF (ABS(FA).LE.ABS(FB)) GO TO 30
S = FB/FA
IF (A.NE.C) GO TO 20
C LINEAR INTERPOLATION
P = (C-B)*S
Q = ONE-S
GO TO 25
C INVERSE QUADRATIC INTERPOLATION
20 Q = FA/FC
R = FB/FC
RONE = R-ONE
P = S*((C-B)*Q*(Q-R)-(B-A)*RONE)
Q = (Q-ONE)*RONE*(S-ONE)
25 IF (P.GT.ZERO) Q = -Q
IF (P.LT.ZERO) P = -P
S = E
E = D
C IF ABS(P/Q).GE.75*ABS(C-B) THEN
C FORCE BISECTION
IF (P+P.GE.THREE*RM*Q) GO TO 30
C IF ABS(P/Q).GE..5*ABS(S) THEN FORCE
C BISECTION. S = THE VALUE OF P/Q
C ON THE STEP BEFORE THE LAST ONE
IF (P+P.GE.ABS(S*Q)) GO TO 30
D = P/Q
GO TO 35
C BISECTION
30 E = RM
D = E
C INCREMENT B
35 A = B
FA = FB
TEMP = D
IF (ABS(TEMP).LE.HALF*TOL) TEMP = SIGN(HALF*TOL,RM)
B = B+TEMP
S = B
FB = F(S)
IC = IC+1
IF (FB*FC.LE.ZERO) GO TO 10
GO TO 5
C CONVERGENCE OF B
40 A = C
MAXFN = IC
GO TO 9005
C MAXFN EVALUATIONS
45 IER = 129
A = C
MAXFN = IC
GO TO 9000
C TERMINAL ERROR - F(A) AND F(B) HAVE
C THE SAME SIGN
50 IER = 130
MAXFN = IC
9000 CONTINUE
c CALL UERTST (IER,6HZBRENT)
9005 RETURN
END
C *******************************************************
C This is an external function subprogram f(x) whose root
C is to be found. This function, in our case, is f(rho) and
C is given by
C f(rho) = Ractual - Rcal
C = Ractual - (R0 + rho_guess*cos(theta_guess
C + delta*sin(theta_guess)))
C *************************************************************
function f(x)
Implicit Real*8 (A-H,O-Z)
Real*8 ellip,Delta,Rmaj,Rmin,ractual,Theta,zactual
logical quadrant
COMMON /COMBL/ rarr,zarr,Rhoo,dRhor,dRhoz
common /combl1/ alpha,omega,vellight,rmaj,rmin,
& peakdensity,zk0,pi
COMMON /COMBL2/ ellip,delta
common /combl3/ ractual,zactual,quadrant,zslope
If (dsqrt(zactual**2+(ractual-rmaj)**2).lt.1.0d-4) then
theta = pi/2
else
Theta = dasin(Zactual/ellip/x)
endif
if (quadrant) theta = pi - theta
f=ractual-(Rmaj+x*dcos(theta+delta*zactual/ellip/x))
return
end
C *************************************************************
C Fourth order Runge-Kutta integrator
C *************************************************************
SUBROUTINE rk4sys(n,t,x,h)
implicit real*8 (a-h,o-z)
dimension x(4),y(4),f1(4),f2(4),f3(4),f4(4)
c print 7,t,(x(i),i=1,n)
h2 = 0.5 * h
CALL fcn(x,f1)
do 2 i = 1,n
y(i) = x(i) + h2*f1(i)
2 continue
CALL fcn(y,f2)
do 3 i = 1,n
y(i) = x(i) + h2*f2(i)
3 continue
CALL fcn(y,f3)
do 4 i = 1,n
y(i) = x(i) + h*f3(i)
4 continue
CALL fcn(y,f4)
do 5 i = 1,n
x(i) = x(i) + h*(f1(i) + 2.0*(f2(i) + f3(i)) + f4(i))/6.0
5 continue
t = t + h
c print 7,t,(x(i),i = 1,n)
c7 format(' RK4. t,x: ',d10.4,2(1x,d10.4))
6 continue
return
END
C *************************************************************
C This Subroutine calculates derivatives of R,Z,k_R & k_z
C *************************************************************
c Subroutine FCN(NEqn,Time, Y,yprime)
Subroutine FCN(Y,yprime)
Parameter(neqn=4)
Implicit Real*8 (A-H,O-Z)
Real*8 Y(NEqn),yprime(NEqn),alpha,omega,vellight,
& rmaj,rmin
common /combl1/ alpha,omega,vellight,rmaj,rmin,
& peakdensity,zk0,pi
Call Gradient(NEqn,Y,dndR,dndz,Density)
Csq = vellight * vellight
yprime(1) = (Csq/Omega) * Y(3) * zk0
yprime(2) = (Csq/Omega) * Y(4) * zk0
yprime(3) = -(Alpha*0.5d0/Omega) * dndR / zk0
yprime(4) = -(Alpha*0.5d0/Omega) * dndz / zk0
c yprime(5) = dndr * yprime(1) + dndz * yprime(2)
return
End
C *************************************************************
C This subroutine calculates n,dn/dR,dn/dz
C *************************************************************
Subroutine Gradient(NEqn,Y,dndR,dndz,Density)
Parameter(npt=85)
Implicit Real*8 (A-h,O-Z)
Real*8 Y(NEqn),rmaj,rmin,zk0,alpha,omega,
& Rat,rho drhodr,drhodz,d1,d2,density,znedg,prexp,dndrho,
& dndr,dndz
Real*8 rarr(npt),zarr(npt),Rhoo(npt,npt),dRhor(npt,npt),
& dRhoz(npt,npt)
common /combl1/ alpha,omega,vellight,rmaj,rmin,
& peakdensity,zk0,pi
COMMON /COMBL/ rarr,zarr,Rhoo,dRhor,dRhoz
common /prexp/ prexp
call Locate(y(1),y(2),rho,drhodr,drhodz)
Znedg = 1.0d18
Rat = rho/rmin
if(rat .gt. 0.999d0) rat = 0.999d0
d1 = 1.0d0 - rat*rat
d2 = -2.0d0*rat/rmin
If (prexp .eq. 0.5d0)then
density = Peakdensity*sqrt(d1) + znedg
If (abs(rat) .ge. 0.999d0)then
dndrho = 0.0d0
else
dndrho=0.5d0*peakdensity*d2
dndrho = dndrho/sqrt(d1)
Endif
Elseif (prexp .eq. 1.0d0)then
density = Peakdensity*d1 + znedg
If (abs(rat) .ge. 0.999d0)then
dndrho = 0.0d0
else
dndrho = peakdensity * d2
Endif
Elseif (prexp .eq. 2.0d0)then
density = Peakdensity*d1 **2+ znedg
If (abs(rat) .ge. 0.999d0)then
dndrho = 0.0d0
else
dndrho=2.0d0*peakdensity*d2*d1
Endif
Endif
dndr = dndrho*drhodr
dndz = dndrho*drhodz
End
C *************************************************************
C This Subroutine calculates the values of rho, drho/dR, drho/dz
C by interpolation between the grid points at the ray location
C *************************************************************
subroutine Locate(r,z,rho,drhodr,drhodz)
parameter(npt=85)
Implicit Real*8 (A-H,O-Z)
Real*8 rarr(npt),zarr(npt),Rhoo(npt,npt),dRhor(npt,npt),
& dRhoz(npt,npt)
Real*8 Rho,dRhodr,dRhodz,TEMP,TEMP1,TEMP2,TEMP3
rEAL*8 dELTAR,DELTAZ,rmaj,rmin,ELLIP,DELTA
Real*8 d1,d2,rhoa,rhob,drhodra,drhodrb,drhodza,drhodzb
COMMON /COMBL/ rarr,zarr,Rhoo,dRhor,dRhoz
common /combl1/ alpha,omega,vellight,rmaj,rmin,
& peakdensity,zk0,pi
COMMON /COMBL2/ ellip,delta
Deltar = rarr(2) - rarr(1)
Deltaz = zarr(2) - zarr(1)
c i1=((R-rmaj+rmin+2*deltar)/deltar)+1
i1=((R-rarr(1))/deltar)+1
I2 = I1+1
c j1=((z+ ellip*rmin+2*deltaz)/deltaz)+1
j1=((z-zarr(1))/deltaz)+1
c Print*,'ellip=',ellip,' rmin=',rmin,' prod=',ellip*rmin
c Print*,'zarr(1)=',zarr(1),' zlow=',ellip*rmin+2*deltaz
c Print*,'zarr(162)=',zarr(162)
c Print*,'zarr(163)=',zarr(163)
c Print*,'zarr(164)=',zarr(164)
J2 = J1 + 1
D1 = (r-rarr(i1))/deltar
D2 = (z - zarr(j1))/deltaz
c Print*,'i1 =',i1
c Print*,'i2 =',i2
c Print*,'j1 =',j1
c Print*,'j2 =',j2
c Print*,'r=',r,' deltar=',deltar
c Print*,'rarr(i1)=',rarr(i1),' rarr(i2)=',rarr(i2)
c Print*,'z=',z,' deltaz=',deltaz
c Print*,'zarr(j1)=',zarr(j1),' zarr(j2)=',zarr(j2)
if(abs(d1) .gt. 1.000001d0 .or. abs(d2) .gt. 1.0d0)then
c print*,'d1,d2=',d1,d2
endif
if(abs(d1) .gt. 1.0d0)then
RARR(I1) = (RARR(I1-1)+RARR(I2))/2
D1 = (r - rarr(I1))/deltaR
ENDIF
if(abs(d2) .gt. 1.0d0)then
ZARR(J1) = (ZARR(J1-1)+ZARR(J2))/2
D2 = (z - zarr(j1))/deltaz
ENDIF
if(abs(d1) .gt. 1.0d0 .or. abs(d2) .gt. 1.0d0)then
if(abs(abs(d1)-1.0d0).lt.1.0d-10)d1=1.0d0
if(abs(abs(d2)-1.0d0).lt.1.0d-10)d2=1.0d0
c print*,'d1,d2=',d1,d2
endif
if(abs(d1) .gt. 1.0d0 .or. abs(d2) .gt. 1.0d0)then
print*,'d1,d2=',d1,d2
stop
endif
Rhoa = Rhoo(i1,j1) + (Rhoo(i1,j2) - Rhoo(i1,j1))*d2
Rhob = Rhoo(i2,j1) + (Rhoo(i2,j2) - Rhoo(i2,j1))*d2
Rho = Rhoa + (Rhob - Rhoa)*d1
dRhodra = dRhor(i1,j1) + (dRhor(i1,j2) - dRhor(i1,j1))*d2
dRhodrb = dRhor(i2,j1) + (dRhor(i2,j2) - dRhor(i2,j1))*d2
dRhodr = dRhodra + (dRhodrb - dRhodra)*d1
dRhodza = dRhoz(i1,j1) + (dRhoz(i1,j2) - dRhoz(i1,j1))*d2
dRhodzb = dRhoz(i2,j1) + (dRhoz(i2,j2) - dRhoz(i2,j1))*d2
dRhodz = dRhodza + (dRhodzb - dRhodza)*d1
END
Write, Run & Share Fortran code online using OneCompiler's Fortran online compiler for free. It's one of the robust, feature-rich online compilers for Fortran language, running on the latest version 7. Getting started with the OneCompiler's Fortran compiler is simple and pretty fast. The editor shows sample boilerplate code when you choose language as Fortran and start coding.
OneCompiler's Fortran online editor supports stdin and users can give inputs to programs using the STDIN textbox under the I/O tab. Following is a sample Fortran program which takes name as input and prints hello message with your name.
program hello
character :: name*30
read *, name
print *, "Hello ", name
end program hello
Fortran language was initially developed for scientific calculations by IBM in 1957. It has a number of in-built functions to perform mathematical calculations and is ideal for applications which has more mathematical calculations.
| Data type | Description | Usage |
|---|---|---|
| Integer | To store integer variables | integer :: x |
| Real | To store float values | real :: x |
| Complex | To store complex numbers | complex :: x,y |
| Logical | To store boolean values True or false | logical :: x=.True. , logical :: x = .FALSE. |
| Character | To store characters and strings | character :: x |
Variable is a name given to the storage area in order to manipulate them in our programs.
data type :: variable_name
Array is a collection of similar data which is stored in continuous memory addresses.
data-type, dimension (x,y) :: array-name
integer, dimension(3,3) :: cube
Do is used to execute a set of statement(s) iteratively when a given condition is true and the loop variable must be an integer.
do i = start, stop [,step]
! code
end do
Do-While is used to execute a set of statement(s) iteratively when a given condition is true.
do while (condition)
!Code
end do
If is used to execute a set of statements based on a condition.
if (logical-expression) then
!Code
end if
If is used to execute a set of statements based on a condition and execute another set of statements present in else block, if condition specified in If block fails.
if (logical-expression) then
!code when the condition is true
else
!code when the condition fails
end if
Case is similar to switch in C language.
[name:] select case (regular-expression)
case (value1)
! code for value 1
... case (value2)
! code for value 2
...
case default
! default code
...
end select [name]