c = inverse(sqrt(epsilon nought *mu nought))
- From: Wade Ward <zaxfuuq@xxxxxxxxx>
- Date: 15 May 2007 15:25:58 -0700
This program is simply not going away:
!********************************************************
! APPLICATION OF THE FINITE DIFFERENCE METHOD
! This program involves the penetration of a lossless dielectric
SPHER
! The program provides in the maximum absolute value of
! Ey and Ez during the final half-wave of time-stepping
! Assumption:
! +y-directed incident wave with components Ez and Hx.
! I,J,K,NN correspond to X,Y,Z, and Time.
! IMAX,JMAX,KMAX are the maximum values of x,y,z
! NNMAX is the total number of timesteps.
! NHW represents one half-wave cycle.
! MED is the number of different uniform media sections.
! JS is the j-position of the plane wave front.
!
! THIS PROGRAM WAS DEVELOPED BY V. BEMMEL [43]
! AND LATER IMPROVED BY D. TERRY
!*******************************************************
PARAMETER( IMAX=19, JMAX=39, KMAX=19)
PARAMETER( NMAX=2, NNMAX=500, NHW=40, MED=2, JS=3 )
PARAMETER( DELTA=3E-3, CL=3.0E8, F=2.5E9 )
PARAMETER( PIE=3.141592654)
integer, parameter :: dp = kind(1.0d0)
! Define scatterer dimensions
PARAMETER( OI=19.5, OJ=20.0, OK=19.0, RADIUS=15.0)
DIMENSION EX(0:IMAX+1, 0:JMAX+1, 0:KMAX+1,
0:NMAX), &
& EY(0:IMAX+1, 0:JMAX+1, 0:KMAX+1,
0:NMAX), &
& EZ(0:IMAX+1, 0:JMAX+1, 0:KMAX+1,
0:NMAX), &
& HX(0:IMAX+1, 0:JMAX+1, 0:KMAX+1,
0:NMAX), &
& HY(0:IMAX+1, 0:JMAX+1, 0:KMAX+1,
0:NMAX), &
& HZ(0:IMAX+1, 0:JMAX+1, 0:KMAX+1,
0:NMAX), &
& EY1(0:JMAX+1), EZ1(0:JMAX
+1), &
& ER(MED), SIG(MED), CA(MED), CB(MED)
DIMENSION IXMED(0:IMAX+1, 0:JMAX+1, 0:KMAX+1)
DIMENSION IYMED(0:IMAX+1, 0:JMAX+1, 0:KMAX+1)
DIMENSION IZMED(0:IMAX+1, 0:JMAX+1, 0:KMAX+1)
! STORE CB(M)/RB HERE
DIMENSION CBMRB( 2 )
! CONSTITUTIVE PARAMETERS
DATA ER/1.0,4.0/
DATA SIG/0.1,0.0/
!
real(kind=dp) :: x
! Statement function to compute position w.r.t. center of the sphere
!
position(RI,RJ,RK)=SQRT((RI-OI)**2+(RJ-OJ)**2+(RK-OK)**2)
E0=(1E-9)/(36*PIE)
U0=(1E-7)*4*PIE
DT=DELTA/(2*CL)
R=DT/E0
RA=(DT**2)/(U0*E0*(DELTA**2))
RB=DT/(U0*DELTA)
TPIFDT = 2.0*PIE*F*DT
x = 0.0_dp
x = (E0*U0)**(-.5)
write (*,*) x
!********************************************************
! EP # 1 - COMPUTE MEDIA PARAMETERS
!********************************************************
DO 1 M=1,MED
CA(M)=1.0-R*SIG(M)/ER(M)
CB(M)=RA/ER(M)
CBMRB(M) = CB(M)/RB
1 CONTINUE
!
! (i) CALCULATE THE REAL/ACTUAL GRID POINTS
!
! Initialize the media arrays.Index (M) determines which
! medium each point is actually located in and is used to
! index into arrays which determine the constitutive
! parameters of the medium.There are separate M determining
! arrays for EX, EY, and EZ. These arrays correlate the
! integer values of I,J,K to the actual position within
! the lattice. Computing these values now and storing them in these
! arrays as opposed to computing them each time they are
! needed saves a large amount of computation time.
!
DO 3 I=0, IMAX+1
DO 3 J=0, JMAX+1
DO 3 K=0, KMAX+1
IF(position(float(I)+0.5,float(J),float(K)).LE.RADIUS) then
IXMED( I, J, K ) = 2
else
IXMED( I, J, K ) = 1
endif
IF(position(float(I),float(J)+0.5,float(K)).LE.RADIUS) then
IYMED( I, J, K ) = 2
else
IYMED( I, J, K ) = 1
endif
IF(position(float(I),float(J),float(K)+0.5).LE.RADIUS) then
IZMED( I, J, K ) = 2
else
IZMED( I, J, K ) = 1
endif
3 CONTINUE
!*************************************************
! STEP # 2 - INITIALIZE FIELD COMPONENTS
!*************************************************
! components for output
DO 4 J=0,JMAX+1
EY1(J)=0
EZ1(J)=0
4 END DO
DO 5 I=0,IMAX+1
DO 5 J=0,JMAX+1
DO 5 K=0,KMAX+1
DO 5 N = 0, NMAX
EX(I,J,K,N)=0
EY(I,J,K,N)=0
EZ(I,J,K,N)=0
HX(I,J,K,N)=0
HY(I,J,K,N)=0
HZ(I,J,K,N)=0
5 CONTINUE
print *, 'initialization complete'
!********************************************************
! STEP # 3 - USE FD/TD ALGORITHM TO GENERATE
! FIELD COMPONENTS
!********************************************************
! SINCE ONLY FIELD COMPONENTS AT CURRENT TIME (t) AND PREVIOUS
! TWO TIME STEPS ( t-1 AND t-2) ARE REQUIRED FOR COMPUTATION,
! WE SAVE MEMORY SPACE BY USING THE FOLLOWING INDICES
! NCUR is index in for time t
! NPR1 is index in for t-1
! NPR2 is index in for t-2
!
NCUR = 2
NPR1 = 1
NPR2 = 0
! TIME LOOP
DO 15 NN = 1, NNMAX
IF (MOD( NN,10).EQ.0) THEN
! DISPLAY PROGRESS
print *,'NN=',NN
ENDIF
! Next time step - move indices up a notch.
NPR2 = NPR1
NPR1 = NCUR
NCUR = MOD( NCUR+1, 3)
! Z LOOP
DO 20 K=0,KMAX
! Y LOOP
DO 25 J=0,JMAX
! X LOOP
DO 30 I=0,IMAX
!
! (ii)-APPLY SOFT LATTICE TRUNCATION CONDITIONS
!
!---x=delta/2
IF(I.EQ.0) THEN
IF ((K.NE.KMAX).AND.(K.NE.0)) THEN
HY(0,J,K,NCUR) =
(HY(1,J,K-1,NPR2) &
& + HY(1,J,K,NPR2)+HY(1,J,K+1,NPR2))/3
HZ(0,J,K,NCUR)=(HZ(1,J,K-1,NPR2) &
& + HZ(1,J,K,NPR2)+HZ(1,J,K+1,NPR2))/3
else
IF (K.EQ.KMAX) THEN
HY(0,J,KMAX,NCUR) =
(HY(1,J,KMAX-1,NPR2) &
& + HY(1,J,KMAX,NPR2))/2.
HZ(0,J,K,NCUR)=( HZ(1,J,K-1,NPR2) &
& + HZ(1,J,K,NPR2) )/2
ELSE
HY(0,J,K,NCUR) =
( HY(1,J,K,NPR2) &
& + HY(1,J,K+1,NPR2))/2.
HZ(0,J,0,NPR2)=(HZ(1,J,
0,NPR2) &
& + HZ(1,J,1,NPR2))/2
ENDIF
ENDIF
ENDIF
!---y=0
IF(J.EQ.0) THEN
EX(I,0,K,NCUR)=EX(I,1,K,NPR2)
EZ(I,0,K,NCUR)=EZ(I,1,K,NPR2)
ELSE
!---y=ymax
IF(J.EQ.JMAX) THEN
EX(I,JMAX,K,NCUR)=EX(I,JMAX-1,K,NPR2)
EZ(I,JMAX,K,NCUR)=EZ(I,JMAX-1,K,NPR2)
ENDIF
ENDIF
!---z=0
IF(K.EQ.0) THEN
IF ((I.NE.0).AND.(I.NE.IMAX)) THEN
EX(I,J,0,NCUR) = (EX(I-1,J,
1,NPR2) &
& + EX(I,J,1,NPR2)+EX(I+1,J,1,NPR2))/3
EY(I,J,0,NCUR) = (EY(I-1,J,
1,NPR2) &
& + EY(I,J,1,NPR2)+EY(I+1,J,1,NPR2))/3
ELSE
IF(I.EQ.0) THEN
EX(0,J,0,NCUR)=(EX(0,J,1,NPR2)+EX(1,J,1,NPR2))/2
EY(I,J,0,NCUR)=(EY(I,J,1,NPR2)+EY(I+1,J,1,NPR2))/2
ELSE
EX(I,J,0,NCUR)=(EX(I-1,J,1,NPR2)+EX(I,J,1,NPR2))/2
EY(I,J,0,NCUR)=(EY(I-1,J,1,NPR2)+EY(I,J,1,NPR2))/2
ENDIF
ENDIF
ENDIF
!
! (iii) APPLY FD/TD ALGORITHM
!
!-----a. HX generation:
HX(I,J,K,NCUR)=HX(I,J,K,NPR1)+RB*(EY(I,J,K
+1,NPR1) &
& - EY(I,J,K,NPR1)+EZ(I,J,K,NPR1)-EZ(I,J+1,K,NPR1))
!-----b. HY generation:
HY(I,J,K,NCUR)=HY(I,J,K,NPR1)+RB*(EZ(I
+1,J,K,NPR1) &
& - EZ(I,J,K,NPR1)+EX(I,J,K,NPR1)-EX(I,J,K+1,NPR1))
!-----c. HZ generation:
HZ(I,J,K,NCUR)=HZ(I,J,K,NPR1)+RB*(EX(I,J
+1,K,NPR1) &
& - EX(I,J,K,NPR1)+EY(I,J,K,NPR1)-EY(I+1,J,K,NPR1))
!---k=kmax ! SYMMETRY
IF(K.EQ.KMAX) THEN
HX(I,J,KMAX,NCUR)=HX(I,J,KMAX-1,NCUR)
HY(I,J,KMAX,NCUR)=HY(I,J,KMAX-1,NCUR)
ENDIF
!-----d. EX generation:
IF((J.NE.0).AND.(J.NE.JMAX).AND.(K.NE.0)) THEN
M = IXMED( I, J, K )
EX(I,J,K,NCUR) = CA(M)*EX(I,J,K,NPR1) +
CBMRB(M)* &
& (HZ(I,J,K,NCUR)-HZ(I,J-1,K,NCUR)
+HY(I,J,K-1,NCUR) &
& -HY(I,J,K,NCUR))
ENDIF
!-----e. EY generation:
IF(K.NE.0) THEN
M = IYMED( I, J, K )
EY(I,J,K,NCUR)=CA(M)*EY(I,J,K,NPR1) +
CBMRB(M)* &
& (HX(I,J,K,NCUR)-HX(I,J,K-1,NCUR)+HZ(I-1,J,K,NCUR) &
& -HZ(I,J,K,NCUR))
! was unpaired right brace^^
ENDIF
!-----f. EZ generation:
IF ((J.NE.0).AND.(J.NE.JMAX)) THEN
M = IZMED( I, J, K )
! sig(ext)=0 for Ez only from Taflove[14]
IF(M.EQ.1) THEN
CAM=1.0
ELSE
CAM=CA(M)
ENDIF
EZ(I,J,K,NCUR)=CAM*EZ(I,J,K,NPR1)+CBMRB(M)* &
& (HY(I,J,K,NCUR)-HY(I-1,J,K,NCUR)
+HX(I,J-1,K,NCUR) &
& -HX(I,J,K,NCUR))
!
! (iv) APPLY THE PLANE-WAVE SOURCE
!
IF(J.EQ.JS) THEN
EZ(I,JS,K,NCUR) = EZ(I,JS,K,NCUR)+SIN( TPIFDT*NN )
ENDIF
ENDIF
!---i=imax+1/2 ! SYMMETRY
IF(I.EQ.IMAX) THEN
EY(IMAX+1,J,K,NCUR)=EY(IMAX,J,K,NCUR)
EZ(IMAX+1,J,K,NCUR)=EZ(IMAX,J,K,NCUR)
ENDIF
!---k=kmax
IF(K.EQ.KMAX) THEN
EX(I,J,KMAX+1,NCUR)=EX(I,J,KMAX-1,NCUR)
EY(I,J,KMAX+1,NCUR)=EY(I,J,KMAX-1,NCUR)
ENDIF
! I LOOP
30 END DO
!********************************************************
! STEP # 4 - RETAIN THE MAXIMUM ABSOLUTE VALUES DURING
! THE LAST HALF-WAVE
!********************************************************
IF ( (K.EQ.KMAX).AND.(NN.GT.(NNMAX-NHW)) ) THEN
TEMP = ABS( EY(IMAX,J,KMAX-1,NCUR) )
IF (TEMP .GT. EY1(J) ) THEN
EY1(J) = TEMP
ENDIF
TEMP = ABS( EZ(IMAX,J,KMAX,NCUR) )
IF (TEMP .GT. EZ1(J) ) THEN
EZ1(J) = TEMP
ENDIF
ENDIF
! J LOOP
25 END DO
! K LOOP
20 END DO
! NN (time) LOOP
15 END DO
!*******************************************************
!-----Output E as a function of j
DO 130 J = 0, JMAX-1
WRITE(6,135) J+1, float(J), EY1(J)
130 END DO
DO 140 J = 0, JMAX-1
WRITE(6,135) J+1, float(J), EZ1(J)
140 END DO
135 FORMAT(I5,2F15.8)
STOP
END
!end source
I wanted to have a clean desktop before I got hooked up to usenet
again and was fiddling with the above source. I saw U0 and E0 in this
source and reflected on what mischief a person could create
therewith. I knew I wanted another double precision real, and
although I'd done it before, I couldn't recall the syntax.
Not to worry. I have old messages on my machine and found:
integer, parameter :: dp = kind(1.0d0)
real(kind=dp) :: x
x = 0.0_dp
This procedure has been recommended by the experts.
Hideous, it looks like fodder for the fools en route to self
destruction. <snip>
! end old message
I added the above, but not without compiler complaint with respect to
*where* the statements go. When I put them before:
position(RI,RJ,RK)=....
, it creates the following error:
: error 381 - POSITION is not an array.
When I place it right after that block, it compiles fine. What gives?
Am I right to think that this gives the speed of light as a function
of pi alone?
--
Wade Ward
.
- Follow-Ups:
- Re: c = inverse(sqrt(epsilon nought *mu nought))
- From: Richard Maine
- Re: c = inverse(sqrt(epsilon nought *mu nought))
- Prev by Date: Re: Fortran to matlab infuriating problem
- Next by Date: Re: pretty output using Fortran
- Previous by thread: passing large arrays to functions
- Next by thread: Re: c = inverse(sqrt(epsilon nought *mu nought))
- Index(es):
Relevant Pages
|
