|
 |
 |
 |
 |
Fortran Programming Language
|
 |
 |
 |
 |
 |
 |
 |
 |
c = inverse(sqrt(epsilon nought *mu nought))
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 ...
read more »
Wade Ward <zaxf@gmail.com> wrote:
> > x = 0.0_dp
...
> 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?
Well, the compiler is right. Position isn't an array. It is intended to
be a statement function, as the comment above it notes.
Statement functions have loist of quirks, which is why many people,
myself included, don't recommend using them in new code. I won't even
try to list all the quirks. It would take way too long (and I'd miss a
bunch of them unless I spent a few hours checking first).
Statement functions are syntactically very easy to confuse with
assignment statements. One result of that is that the compiler can get
confused also. If you have something that isn't a valid statement
function, you can get an error message from the compiler trying to
interpret it as an assignment. That's what happened here. The compiler
guessed that this might be intended as an assignment to an element of an
array named position, but there is no such array declared, so you get an
error message.
It isn't a valid statement function when you put the x=0.0_dp before it
because x=0.0_dp is an executable statement. Statement functions must be
in the specification part - i.e. before any executable statements.
What puzzles me is why you even have x in this code. It seems to have
nothing to do with the rest of the code; nothing ever references it.
Since it doesn't have anything to do with the rest of the code, it is a
bit hard to suggest sensibly what to do with it.
--
Richard Maine | Good judgement comes from experience;
email: last name at domain . net | experience comes from bad judgement.
domain: summertriangle | -- Mark Twain
On May 15, 6:53 pm, nos@see.signature (Richard Maine) wrote:
> Wade Ward <zaxf@gmail.com> wrote:
> > > x = 0.0_dp
> ...
> > 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?
> Well, the compiler is right. Position isn't an array. It is intended to
> be a statement function, as the comment above it notes.
> Statement functions have loist of quirks, which is why many people,
> myself included, don't recommend using them in new code. I won't even
> try to list all the quirks. It would take way too long (and I'd miss a
> bunch of them unless I spent a few hours checking first).
> Statement functions are syntactically very easy to confuse with
> assignment statements. One result of that is that the compiler can get
> confused also. If you have something that isn't a valid statement
> function, you can get an error message from the compiler trying to
> interpret it as an assignment. That's what happened here. The compiler
> guessed that this might be intended as an assignment to an element of an
> array named position, but there is no such array declared, so you get an
> error message.
> It isn't a valid statement function when you put the x=0.0_dp before it
> because x=0.0_dp is an executable statement. Statement functions must be
> in the specification part - i.e. before any executable statements.
> What puzzles me is why you even have x in this code. It seems to have
> nothing to do with the rest of the code; nothing ever references it.
> Since it doesn't have anything to do with the rest of the code, it is a
> bit hard to suggest sensibly what to do with it.
Statement functions? That's a new one. I'll have a look at MR&C.
x was set to the LHS of statement that involved U0 and E0. Since it
was buried in 6k of stuff that I didn't know how to trim out, it's
understandable that a person would miss the statement, which was:
x = (E0*U0)**(-.5)
EO was 10**-7 * 4pi and UO was 10**-9/36 pi . In the product, the
pi's cancel leaving the inverse square root of (10**-16 * 1/9). Both
are perfect squares, ergo
x= 1/(10**-8 * 1/3) which leaves x at 3 x 10**8. The speed of light
is as much of a constant as a person can imagine.
--
Wade Ward
Wade Ward <zaxf@gmail.com> wrote:
> Statement functions? That's a new one. I'll have a look at MR&C.
More like an old one. If you are looking in MR&C, be sure to check their
section on old features; it is probably buried in there.
--
Richard Maine | Good judgement comes from experience;
email: last name at domain . net | experience comes from bad judgement.
domain: summertriangle | -- Mark Twain
On May 17, 8:39 pm, Wade Ward <zaxf@gmail.com> wrote:
> On May 15, 6:53 pm, nos@see.signature (Richard Maine) wrote:
> > Wade Ward <zaxf@gmail.com> wrote:
> > > > x = 0.0_dp
> > ...
> > > 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?
> > Well, the compiler is right. Position isn't an array. It is intended to
> > be a statement function, as the comment above it notes.
> > Statement functions have loist of quirks, which is why many people,
> > myself included, don't recommend using them in new code. I won't even
> > try to list all the quirks. It would take way too long (and I'd miss a
> > bunch of them unless I spent a few hours checking first).
> > Statement functions are syntactically very easy to confuse with
> > assignment statements. One result of that is that the compiler can get
> > confused also. If you have something that isn't a valid statement
> > function, you can get an error message from the compiler trying to
> > interpret it as an assignment. That's what happened here. The compiler
> > guessed that this might be intended as an assignment to an element of an
> > array named position, but there is no such array declared, so you get an
> > error message.
> > It isn't a valid statement function when you put the x=0.0_dp before it
> > because x=0.0_dp is an executable statement. Statement functions must be
> > in the specification part - i.e. before any executable statements.
> > What puzzles me is why you even have x in this code. It seems to have
> > nothing to do with the rest of the code; nothing ever references it.
> > Since it doesn't have anything to do with the rest of the code, it is a
> > bit hard to suggest sensibly what to do with it.
> Statement functions? That's a new one. I'll have a look at MR&C.
> x was set to the LHS of statement that involved U0 and E0. Since it
> was buried in 6k of stuff that I didn't know how to trim out, it's
> understandable that a person would miss the statement, which was:
> x = (E0*U0)**(-.5)
> EO was 10**-7 * 4pi and UO was 10**-9/36 pi . In the product, the
> pi's cancel leaving the inverse square root of (10**-16 * 1/9). Both
> are perfect squares, ergo
> x= 1/(10**-8 * 1/3) which leaves x at 3 x 10**8. The speed of light
> is as much of a constant as a person can imagine.
> --
> Wade Ward
Perhaps it might help to look at the text where these examples may
have originated. Based on the exercise numbers of two programs cited
and chapter titles (finite difference and method of moments) along
with Appendix D (it covers solution of simultaneous equations), I
suspect these may have come from "Numerical Techniques in
Electromagnetics", Matthew N. O. Sadiku, 2000, CRC Press, ISBN
0849313953. [Thanks to Google book search.]
This still leaves un-resolved this program's problems with out-of-
bounds references to arrays. In addition, there is quite a bit of
archaic code.
In the end, I believe that you need some knowledge of the problem
domain to get a handle on this program. (Sorry, not my field at
all. :-))
-- elliot
"Wade Ward" <zaxf@gmail.com> wrote in message
news:1179448796.877248.278660@u30g2000hsc.googlegroups.com...
> x= 1/(10**-8 * 1/3) which leaves x at 3 x 10**8. The speed of light
> is as much of a constant as a person can imagine.
> --
Indeed. It is defined to be 2.99792458e8 m/s exactly.
Regards,
Mike Metcalf
On May 17, 8:39 pm, Wade Ward <zaxf@gmail.com> wrote:
> On May 15, 6:53 pm, nos@see.signature (Richard Maine) wrote:
> > Wade Ward <zaxf@gmail.com> wrote:
> > > > x = 0.0_dp
> > ...
> > > 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?
> > Well, the compiler is right. Position isn't an array. It is intended to
> > be a statement function, as the comment above it notes.
> > Statement functions have loist of quirks, which is why many people,
> > myself included, don't recommend using them in new code. I won't even
> > try to list all the quirks. It would take way too long (and I'd miss a
> > bunch of them unless I spent a few hours checking first).
> > Statement functions are syntactically very easy to confuse with
> > assignment statements. One result of that is that the compiler can get
> > confused also. If you have something that isn't a valid statement
> > function, you can get an error message from the compiler trying to
> > interpret it as an assignment. That's what happened here. The compiler
> > guessed that this might be intended as an assignment to an element of an
> > array named position, but there is no such array declared, so you get an
> > error message.
> > It isn't a valid statement function when you put the x=0.0_dp before it
> > because x=0.0_dp is an executable statement. Statement functions must be
> > in the specification part - i.e. before any executable statements.
> > What puzzles me is why you even have x in this code. It seems to have
> > nothing to do with the rest of the code; nothing ever references it.
> > Since it doesn't have anything to do with the rest of the code, it is a
> > bit hard to suggest sensibly what to do with it.
> Statement functions? That's a new one. I'll have a look at MR&C.
> x was set to the LHS of statement that involved U0 and E0. Since it
> was buried in 6k of stuff that I didn't know how to trim out, it's
> understandable that a person would miss the statement, which was:
> x = (E0*U0)**(-.5)
> EO was 10**-7 * 4pi and UO was 10**-9/36 pi . In the product, the
> pi's cancel leaving the inverse square root of (10**-16 * 1/9). Both
> are perfect squares, ergo
> x= 1/(10**-8 * 1/3) which leaves x at 3 x 10**8. The speed of light
> is as much of a constant as a person can imagine.
> --
> Wade Ward
Indeed, by doing some algebra, some relevant constants take on
interesting values:
[let c = Cl = speed of light]
Eo Uo = 1 / c^2
c Uo = 1 / c Eo
Dt / delta = 1 / 2c
R = delta c Uo /2
RA/RB = c Uo /2
R = 180 pi / 1000
RA = 1/4
RA/RB = 60 pi
Constants used later in the program take on the values:
RB = 1 / (240 pi)
CA = [ 1 - (180 pi)/10000 , 1 - 0 ]
CBMRB = [ 60 pi, 15 pi ]
TPIFDT = (2 pi)/80
Adding a "guard" plane for I=-1 in HY and HZ and setting this to zero
makes very little difference in the leading digits of the output. This
may not be a legitimate way to deal with the out of bounds program but
it is tempting to think that it does not perturb the solution very
much.
-- e-mail: epc8 at juno dot com
On May 18, 7:09 am, e p chandler <e@juno.com> wrote:
> On May 17, 8:39 pm, Wade Ward <zaxf@gmail.com> wrote:
> > On May 15, 6:53 pm, nos@see.signature (Richard Maine) wrote:
> > > Wade Ward <zaxf@gmail.com> wrote:
> > > > > x = 0.0_dp
> > > ...
> > > > 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?
> > > Well, the compiler is right. Position isn't an array. It is intended to
> > > be a statement function, as the comment above it notes.
> > > Statement functions have loist of quirks, which is why many people,
> > > myself included, don't recommend using them in new code. I won't even
> > > try to list all the quirks. It would take way too long (and I'd miss a
> > > bunch of them unless I spent a few hours checking first).
> > > Statement functions are syntactically very easy to confuse with
> > > assignment statements. One result of that is that the compiler can get
> > > confused also. If you have something that isn't a valid statement
> > > function, you can get an error message from the compiler trying to
> > > interpret it as an assignment. That's what happened here. The compiler
> > > guessed that this might be intended as an assignment to an element of an
> > > array named position, but there is no such array declared, so you get an
> > > error message.
> > > It isn't a valid statement function when you put the x=0.0_dp before it
> > > because x=0.0_dp is an executable statement. Statement functions must be
> > > in the specification part - i.e. before any executable statements.
> > > What puzzles me is why you even have x in this code. It seems to have
> > > nothing to do with the rest of the code; nothing ever references it.
> > > Since it doesn't have anything to do with the rest of the code, it is a
> > > bit hard to suggest sensibly what to do with it.
> > Statement functions? That's a new one. I'll have a look at MR&C.
> > x was set to the LHS of statement that involved U0 and E0. Since it
> > was buried in 6k of stuff that I didn't know how to trim out, it's
> > understandable that a person would miss the statement, which was:
> > x = (E0*U0)**(-.5)
> > EO was 10**-7 * 4pi and UO was 10**-9/36 pi . In the product, the
> > pi's cancel leaving the inverse square root of (10**-16 * 1/9). Both
> > are perfect squares, ergo
> > x= 1/(10**-8 * 1/3) which leaves x at 3 x 10**8. The speed of light
> > is as much of a constant as a person can imagine.
> > --
> > Wade Ward
> Perhaps it might help to look at the text where these examples may
> have originated. Based on the exercise numbers of two programs cited
> and chapter titles (finite difference and method of moments) along
> with Appendix D (it covers solution of simultaneous equations), I
> suspect these may have come from "Numerical Techniques in
> Electromagnetics", Matthew N. O. Sadiku, 2000, CRC Press, ISBN
> 0849313953. [Thanks to Google book search.]
> This still leaves un-resolved this program's problems with out-of-
> bounds references to arrays. In addition, there is quite a bit of
> archaic code.
> In the end, I believe that you need some knowledge of the problem
> domain to get a handle on this program. (Sorry, not my field at
> all. :-))
> -- elliot- Hide quoted text -
>
program elliot6
PARAMETER( PIE=3.141592654)
integer, parameter :: dp = kind(1.0d0)
!contains real function position
! Define scatterer dimensions
PARAMETER( OI=19.5, OJ=20.0, OK=19.0, RADIUS=15.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)
! another statement function
asf(arg1) = SQRT(arg1-OI)**2
E0=(1E-9)/(36*PIE)
U0=(1E-7)*4*PIE
C =(E0*U0)**(-.5)
x = 0.0_dp
x = (E0*U0)**(-.5)
write (*,*) x, C, E0, U0
i = 4
j = 5
k = 6
IF(position(float(I)+0.5,float(J),float(K)).LE.20) then
write (*,*) "tja"
else
write (*,*) "hja"
endif
!$$$$$$ IF( asf(float(I)).GE.20) then
!$$$$$$ write (*,*) "tja"
!$$$$$$ else
!$$$$$$ write (*,*) "hja"
!$$$$$$ endif
end program elliot6
I intend to post this using the google portal and intend it as a
response to Richard Maine and Mike Metcalf as well. Google doesn't
let a person see a tree of messages that one sees in a newsbrowser, so
it's harder to tailor responses and keep track of subthreads.
I found the option in silverfrost to "Integrate with SDBG," yet don't
find a difference in the debugger as I experience it. I have a
variable list but no idea what type of function position or asf is. I
can't step into them. So I'm a little stuck guessing how to replace
them with an explicit function.
I tried to make a one-dimensional analog of position with asf. What I
have compiles but causes my computer to beep on execution instead of
giving the other branch of the if-then-else control. The variables
that serve as the arguments for the statement function remain
undefined for the duration of the program, that is, RI,RJ,RK and
arg1. Should I expect otherwise?
I've always lobbied for the existence of C-style comments in fortran,
but found that my IDE had a way to comment out (and in) highlighted
blocks, as seen in the final if statement. Very handy stuff. Right
next to those buttons are ones that will decrease (increase) the
spacing for an entire block.
With MM posting a value for c at tens of thousands of meters per
second less than 3 x 10**8, this forces either
C =(E0*U0)**(-.5) to be an approximation or
E0=(1E-9)/(36*PIE)
U0=(1E-7)*4*PIE are approximations. While I took undergrad
electromagnetics, we never calculated epsilon nought or mu nought.
These were constants in Maxwell's equations that no one argued about.
My guess is that to solve for them, you would need to make a contrived
choice of parameters and gaussian surfaces that isolates them during
the integration. The formula for the inverse of the square root of
their product equalling the speed of light I got from a book I was
reading a few years back by Bryan Green. It may have been _The
Elegant Universe_. Since the formula was in the footnotes, which was
reserved for more technical treatment, I have to think it was right.
This inclines me to the view that the transmittivity and permittivity
in free space are off. I've got a guess here though. In the fourth
of Maxwell's equation there was a term that Dr. Clerk Maxwell added to
make the theory robust. I think this term was ignored when E0 and U0
were calculated. As I recall, it involved the product of E0 and U0
and would have been orders of magnitude smaller, so ignoring it for
the sake of getting a quick and dirty result is understandable.
Anyways, I'll see if a library in Western Michigan has the book you
mention.
--
Wade Ward
In article <TAi3i.18$ix.2@trndny01>,
Michael Metcalf <michaelmetc@compuserve.com> wrote:
>"Wade Ward" <zaxf@gmail.com> wrote in message
>> The speed of light is as much of a constant as a person can imagine.
>Indeed. It is defined to be 2.99792458e8 m/s exactly.
You won't get that many figures with most processors. If you need
to do physical theory with high accuracy 2.99792458d8 is safer.
-- John Harper, School of Mathematics, Statistics and Computer Science,
Victoria University, PO Box 600, Wellington 6140, New Zealand
e-mail john.har@vuw.ac.nz phone (+64)(4)463 5341 fax (+64)(4)463 5045
In article <1179703684.134338.96@n15g2000prd.googlegroups.com>,
Wade Ward <zaxf@gmail.com> wrote:
>With MM posting a value for c at tens of thousands of meters per
>second less than 3 x 10**8, this forces either
>C =(E0*U0)**(-.5) to be an approximation or
>E0=(1E-9)/(36*PIE)
>U0=(1E-7)*4*PIE are approximations.
Not quite. You don't say what units you are using, but in SI units
both C and U0 above are both exact, E0 is an approximation.
The SI definitions give the units of U0 as N/A**2 = kg m s**(-2) A**(-2),
C as m/s, E0 as F/m = kg**(-1) m**(-3) s**4 A**2.
(If you invent suitable derived types you can even make Fortran do
your units as well as your numbers)
See, for example, http://www.physicstoday.org/guide/fundcon.html
-- John Harper, School of Mathematics, Statistics and Computer Science,
Victoria University, PO Box 600, Wellington 6140, New Zealand
e-mail john.har@vuw.ac.nz phone (+64)(4)463 5341 fax (+64)(4)463 5045
John Harper wrote:
> In article <1179703684.134338.96@n15g2000prd.googlegroups.com>,
> Wade Ward <zaxf@gmail.com> wrote:
>> With MM posting a value for c at tens of thousands of meters per
>> second less than 3 x 10**8, this forces either
>> C =(E0*U0)**(-.5) to be an approximation or
>> E0=(1E-9)/(36*PIE)
>> U0=(1E-7)*4*PIE are approximations.
> Not quite. You don't say what units you are using, but in SI units
> both C and U0 above are both exact, E0 is an approximation.
> The SI definitions give the units of U0 as N/A**2 = kg m s**(-2) A**(-2),
> C as m/s, E0 as F/m = kg**(-1) m**(-3) s**4 A**2.
> (If you invent suitable derived types you can even make Fortran do
> your units as well as your numbers)
> See, for example, http://www.physicstoday.org/guide/fundcon.html
Or for a selectable list (with equations where appropriate):
http://physics.nist.gov/cuu/Constants/index.html
cheers,
paulv
--
Paul van Delst Ride lots.
CIMSS @ NOAA/NCEP/EMC Eddy Merckx
On May 20, 9:32 pm, har@mcs.vuw.ac.nz (John Harper) wrote:
> In article <1179703684.134338.96@n15g2000prd.googlegroups.com>,
> Wade Ward <zaxf@gmail.com> wrote:
> >With MM posting a value for c at tens of thousands of meters per
> >second less than 3 x 10**8, this forces either
> >C =(E0*U0)**(-.5) to be an approximation or
> >E0=(1E-9)/(36*PIE)
> >U0=(1E-7)*4*PIE are approximations.
> Not quite. You don't say what units you are using, but in SI units
> both C and U0 above are both exact, E0 is an approximation.
> The SI definitions give the units of U0 as N/A**2 = kg m s**(-2) A**(-2),
> C as m/s, E0 as F/m = kg**(-1) m**(-3) s**4 A**2.
> (If you invent suitable derived types you can even make Fortran do
> your units as well as your numbers)
> See, for example,http://www.physicstoday.org/guide/fundcon.html
I'm able for the first time in a while to surf the net and see at the
above link that c is defined as a rational slightly less than 3 x
10**8 and U0 is indeed as above. Wouldn't that force E0 to be
irrational as well?
--
Wade Ward
Wade Ward wrote:
> On May 20, 9:32 pm, har@mcs.vuw.ac.nz (John Harper) wrote:
>> In article <1179703684.134338.96@n15g2000prd.googlegroups.com>,
>> Wade Ward <zaxf@gmail.com> wrote:
>>> With MM posting a value for c at tens of thousands of meters per
>>> second less than 3 x 10**8, this forces either
>>> C =(E0*U0)**(-.5) to be an approximation or
>>> E0=(1E-9)/(36*PIE)
>>> U0=(1E-7)*4*PIE are approximations.
>> Not quite. You don't say what units you are using, but in SI units
>> both C and U0 above are both exact, E0 is an approximation.
>> The SI definitions give the units of U0 as N/A**2 = kg m s**(-2) A**(-2),
>> C as m/s, E0 as F/m = kg**(-1) m**(-3) s**4 A**2.
>> (If you invent suitable derived types you can even make Fortran do
>> your units as well as your numbers)
>> See, for example,http://www.physicstoday.org/guide/fundcon.html
> I'm able for the first time in a while to surf the net and see at the
> above link that c is defined as a rational slightly less than 3 x
> 10**8 and U0 is indeed as above. Wouldn't that force E0 to be
> irrational as well?
Why is the irrationality a problem? They're still constants:
c = 299792458 m.s^-1 (exact)
u0 = 4pi x 10^-7 N.A^-2 (exact)
e0 = 1/(u0.c^2) F.m^-1 (exact)
where the "exact" refers to the uncertainty of the value (in a vacuum). Since they are all
defined quantities, not measured, their knowledge is considered "exact".
cheers,
paulv
--
Paul van Delst Ride lots.
CIMSS @ NOAA/NCEP/EMC Eddy Merckx
In article <f2sr4o$p8@news.nems.noaa.gov>,
Paul van Delst <Paul.vanDe@noaa.gov> wrote:
>Wade Ward wrote:
>> On May 20, 9:32 pm, har@mcs.vuw.ac.nz (John Harper) wrote:
>>> In article <1179703684.134338.96@n15g2000prd.googlegroups.com>,
>>> Wade Ward <zaxf@gmail.com> wrote:
>>>> C =(E0*U0)**(-.5) to be an approximation or
>>>> E0=(1E-9)/(36*PIE)
>>>> U0=(1E-7)*4*PIE are approximations.
>>> Not quite. You don't say what units you are using, but in SI units
>>> both C and U0 above are both exact, E0 is an approximation.
>>> See, for example,http://www.physicstoday.org/guide/fundcon.html
>> I'm able for the first time in a while to surf the net and see at the
>> above link that c is defined as a rational slightly less than 3 x
>> 10**8 and U0 is indeed as above. Wouldn't that force E0 to be
>> irrational as well?
>Why is the irrationality a problem? They're still constants:
>c = 299792458 m.s^-1 (exact)
>u0 = 4pi x 10^-7 N.A^-2 (exact)
>e0 = 1/(u0.c^2) F.m^-1 (exact)
>where the "exact" refers to the uncertainty of the value (in a vacuum).
>Since they are all defined quantities, not measured, their knowledge
>is considered "exact".
I have another problem with Wade's terminology. c is a velocity, which
happens to be a rational (for that matter an integer) number of metres
per second, which Paul gave. c is NOT a rational number. Rational
numbers are dimensionless. Remember the spacecraft that crashed on Mars
because someone was careless with units...
-- John Harper, School of Mathematics, Statistics and Computer Science,
Victoria University, PO Box 600, Wellington 6140, New Zealand
e-mail john.har@vuw.ac.nz phone (+64)(4)463 5341 fax (+64)(4)463 5045
John Harper wrote:
(snip)
> I have another problem with Wade's terminology. c is a velocity, which
> happens to be a rational (for that matter an integer) number of metres
> per second, which Paul gave. c is NOT a rational number. Rational
> numbers are dimensionless. Remember the spacecraft that crashed on Mars
> because someone was careless with units...
I would agree with you, except that this is the Fortran newsgroup
where variables do not have units. The quantity stored in the
variable is the desired quantity divided by the appropriate units.
It seems to me that physicists and related scientists tend to
consider variables as holding a quantity with units, while
engineers tend to consider the units as part of the equation.
For the latter case, you might see in a book:
F(Newtons) = m (kilograms) * a (meters/second**2).
-- glen
> I'm able for the first time in a while to surf the net and see at the
> above link that c is defined as a rational slightly less than 3 x
> 10**8 and U0 is indeed as above.
IIRC, u0 is just from geometry - the surface area of a unit sphere in some
suitable space. As to c, among it, the meter and the second, you can define a
value for two of them, and the third results from a measurement based on those
standard values. It used to be the case that the meter and the second were
given definitions, and you needed to measure c. But the definition of the
meter is quite inexact: Being based on the wavelength of some line in a
spectrum, its relative error was fairly high, compared to the second (I
remember three or four orders of magnitude). So what was done was to retain
the definition of the second, _define_ c to be a certain value (derived from
the best _measured_ value of c before this decision), and to derive the value
of one meter from those two. The meter is now defined to the same accuracy as
the second, and c is now an exact fundamental constant (as it should be, IMO).
Jan
On May 21, 7:49 pm, glen herrmannsfeldt <g@ugcs.caltech.edu> wrote:
> John Harper wrote:
> (snip)
> > I have another problem with Wade's terminology. c is a velocity, which
> > happens to be a rational (for that matter an integer) number of metres
> > per second, which Paul gave. c is NOT a rational number. Rational
> > numbers are dimensionless. Remember the spacecraft that crashed on Mars
> > because someone was careless with units...
> I would agree with you, except that this is the Fortran newsgroup
> where variables do not have units. The quantity stored in the
> variable is the desired quantity divided by the appropriate units.
> It seems to me that physicists and related scientists tend to
> consider variables as holding a quantity with units, while
> engineers tend to consider the units as part of the equation.
> For the latter case, you might see in a book:
> F(Newtons) = m (kilograms) * a (meters/second**2).
program elliot7
implicit none
integer, parameter :: dp = kind(1.0d0)
real(kind=dp):: x,U0,E0,PIE, c, e0_calc
PIE = 4* atan(1.0)
U0=(1E-7_dp)*4.0 * PIE
E0= 8.854187E-12_dp
c = 2.99792458E8_dp
x = (E0*U0)**(-.5)
write (*,'(f30.20)') x, E0, U0, PIE
write (*, *) "going the other way"
e0_calc = 1/ (U0 * c**2.0)
write (*,'(f35.25)') c, U0, PIE, e0_calc
write (*,*) "12345678.0000000000012345678911234"
write (*,*) "12345678.1234567891123456789212345"
! end source begin output
end program elliot7
299792467.67059397698000000000
0.00000000000885418700
0.00000125663709640503
3.14159274101257324230
going the other way
299792458.0000000000000000000000000
0.0000012566370964050291944
3.1415927410125732423000000
0.0000000000088541875712302
12345678.0000000000012345678911234
12345678.1234567891123456789212345
This is the latest version of this program. I tried to get as much
width as I could for e0_calc. When I took E0 from a website that gave
it to 7 figures it gives the speed of light off by 9 m/s, which is
faster than Ben Johnson. I'd be curious to know how much accuracy I
got.
One question: If I have a variable declared as dp and assign it a
value like (1E-7)*4.0 , am I well-advised to append _dp on the things
on the RHS?
--
Wade Ward
On May 21, 3:15 pm, Paul van Delst <Paul.vanDe@noaa.gov> wrote:
> Wade Ward wrote:
> > On May 20, 9:32 pm, har@mcs.vuw.ac.nz (John Harper) wrote:
> >> In article <1179703684.134338.96@n15g2000prd.googlegroups.com>,
> >> Wade Ward <zaxf@gmail.com> wrote:
> >>> With MM posting a value for c at tens of thousands of meters per
> >>> second less than 3 x 10**8, this forces either
> >>> C =(E0*U0)**(-.5) to be an approximation or
> >>> E0=(1E-9)/(36*PIE)
> >>> U0=(1E-7)*4*PIE are approximations.
> >> Not quite. You don't say what units you are using, but in SI units
> >> both C and U0 above are both exact, E0 is an approximation.
> >> The SI definitions give the units of U0 as N/A**2 = kg m s**(-2) A**(-2),
> >> C as m/s, E0 as F/m = kg**(-1) m**(-3) s**4 A**2.
> >> (If you invent suitable derived types you can even make Fortran do
> >> your units as well as your numbers)
> >> See, for example,http://www.physicstoday.org/guide/fundcon.html
> > I'm able for the first time in a while to surf the net and see at the
> > above link that c is defined as a rational slightly less than 3 x
> > 10**8 and U0 is indeed as above. Wouldn't that force E0 to be
> > irrational as well?
> Why is the irrationality a problem? They're still constants:
> c = 299792458 m.s^-1 (exact)
> u0 = 4pi x 10^-7 N.A^-2 (exact)
> e0 = 1/(u0.c^2) F.m^-1 (exact)
> where the "exact" refers to the uncertainty of the value (in a vacuum). Since they are all
> defined quantities, not measured, their knowledge is considered "exact".
The irrationality of mu nought and epsilon nought isn't a problem, it
simply is. I disagree with John that a number loses rationality or
lack thereof if it has units.
I tried to find a treatment of how u0 is derived. (I miss Halliday
and Resnick.) It doesn't help that it's pronuounced mu nought and
looks like u zero. I think Jan has the right idea, thinking that a
person looks at the unit sphere somehow, but I'm way to rusty to work
it up from first principles.
--
Wade Ward
Wade Ward wrote:
(snip)
> I tried to find a treatment of how u0 is derived. (I miss Halliday
> and Resnick.) It doesn't help that it's pronuounced mu nought and
> looks like u zero. I think Jan has the right idea, thinking that a
> person looks at the unit sphere somehow, but I'm way to rusty to work
> it up from first principles.
Electromagnetic units are derived either through the magnetic
effect (force per unit length between two wires carrying
a specified current a specified distance apart) or electrostatic
(force between two charge distributions a specified distance apart).
MKSA (the A being Ampere) is magnetic based, such that mu nought
is based on geometry and to make the units have a reasonable size.
That is, the Ampere is based on the force between wires with
a specific geometry which sets the value of mu nought in the
appropriate equations.
Epsilon nought is then determined through Maxwell's equations
to satisfy the previously indicated relation.
The more commonly used CGS units are based on setting the
electrical charge unit to specify the force between two charges,
which again results in 4pi in the equations. Another set of
units, Heaviside-Lorentz units, are similar but with the 4pi
moved around such that it doesn't appear in most equations.
But yes, the 4pi comes from a surface integral of the electric
field divided by radius squared, and the area of a sphere
being 4 pi r**2.
-- glen
In article <1179875920.294189.65@w5g2000hsg.googlegroups.com>,
Wade Ward <zaxf@gmail.com> wrote:
>integer, parameter :: dp = kind(1.0d0)
>real(kind=dp):: x,U0,E0,PIE, c, e0_calc
>PIE = 4* atan(1.0)
>U0=(1E-7_dp)*4.0 * PIE
>E0= 8.854187E-12_dp
>c = 2.99792458E8_dp
and then lots of calculations showing not very many correct figures.
With only 7 significant figures quoted in E0 you can't hope for better
in anything depending on it, and atan(1.0) is the single-precision
version of atan. So I'd leave c as is, but change the statements
evaluating PIE, U0 and E0 to
PIE = 4.0_dp*atan(1.0_dp)
U0 = 4e-7_dp * PIE
E0 = 1/ (U0*c**2) ! because **2 is better for integer powers
>One question: If I have a variable declared as dp and assign it a
>value like (1E-7)*4.0 , am I well-advised to append _dp on the things
>on the RHS?
Yes. Two reasons.
1. You sometimes lose accuracy if you don't, so it's safer to get
into the habit of avoiding single-precision real (or complex) numbers
in double-precision expressions.
2. If you think that single precision might be enough or that quadruple
precision is needed (and available), you have to change only one
statement in the entire program: the definition of dp. That's also a
reason to avoid writing PIE = 4d0*atan(1d0) even though it's briefer.
-- John Harper, School of Mathematics, Statistics and Computer Science,
Victoria University, PO Box 600, Wellington 6140, New Zealand
e-mail john.har@vuw.ac.nz phone (+64)(4)463 5341 fax (+64)(4)463 5045
> This is the latest version of this program. I tried to get as much
> width as I could for e0_calc. When I took E0 from a website that gave
> it to 7 figures it gives the speed of light off by 9 m/s, which is
> faster than Ben Johnson. I'd be curious to know how much accuracy I
> got.
Ben Johnson et al. are definitely faster than 9 m/s - they do the 100 m in
less than ten seconds!
As u0 and c are now exact, so should e0 in the same system of units. And the
NIST web site says so. It just doesn't give enough significant figures for
your use.
Jan
On May 23, 2:09 am, Jan Vorbrggen <jvorbrueg@not-mediasec.de>
wrote:
> > This is the latest version of this program. I tried to get as much
> > width as I could for e0_calc. When I took E0 from a website that gave
> > it to 7 figures it gives the speed of light off by 9 m/s, which is
> > faster than Ben Johnson. I'd be curious to know how much accuracy I
> > got.
> Ben Johnson et al. are definitely faster than 9 m/s - they do the 100 m in
> less than ten seconds!
I realized that right after I posted. I named the only people on the
planet who _can_ exceed 9 m/s without machines or free fall or
catastrophic consequences.
--
Wade Ward
Wade Ward <zaxf@gmail.com> wrote:
> On May 23, 2:09 am, Jan Vorbrggen <jvorbrueg@not-mediasec.de>
> wrote:
> > > it gives the speed of light off by 9 m/s, which is
> > > faster than Ben Johnson.
> > Ben Johnson et al. are definitely faster than 9 m/s - they do the 100 m in
> > less than ten seconds!
> I realized that right after I posted. I named the only people on the
> planet who _can_ exceed 9 m/s without machines or free fall or
> catastrophic consequences.
Well, it all depends on your frame of reference. Since the speed of
light was mentioned, seems like relativity ought to be important. I
figure I'm doing something on the order of 300 m/s sitting here at my
desk rotating around the center of the earth. Or quite a bit more around
the Sun or the center of the galaxy. It's enough to make me giddy, or is
that just from my age? :-)
--
Richard Maine | Good judgement comes from experience;
email: last name at domain . net | experience comes from bad judgement.
domain: summertriangle | -- Mark Twain
On May 22, 10:51 pm, har@mcs.vuw.ac.nz (John Harper) wrote:
> and then lots of calculations showing not very many correct figures.
> With only 7 significant figures quoted in E0 you can't hope for better
> in anything depending on it, and atan(1.0) is the single-precision
> version of atan. So I'd leave c as is, but change the statements
> evaluating PIE, U0 and E0 to
> PIE = 4.0_dp*atan(1.0_dp)
> U0 = 4e-7_dp * PIE
> E0 = 1/ (U0*c**2) ! because **2 is better for integer powers
Does an integer exponent make for a mixed-mode calculation?
> >One question: If I have a variable declared as dp and assign it a
> >value like (1E-7)*4.0 , am I well-advised to append _dp on the things
> >on the RHS?
> Yes. Two reasons.
> 1. You sometimes lose accuracy if you don't, so it's safer to get
> into the habit of avoiding single-precision real (or complex) numbers
> in double-precision expressions.
> 2. If you think that single precision might be enough or that quadruple
> precision is needed (and available), you have to change only one
> statement in the entire program: the definition of dp. That's also a
> reason to avoid writing PIE = 4d0*atan(1d0) even though it's briefer.
program elliot8
implicit none
integer, parameter :: dp = kind(1.0d0)
real(kind=dp):: U0,E0,PIE, c, e0_calc
PIE = 4.0_dp* atan(1.0_dp)
U0=(1E-7_dp)*4.0_dp * PIE
E0= 8.854187E-12_dp
c = 2.99792458E8_dp
e0_calc = 1.0_dp/ (U0 * c**2.0_dp)
write (*,'(f35.25)') c, U0, PIE,E0, e0_calc
write (*,*) "12345678.0000000000012345678911234"
write (*,*) "12345678.1234567891123456789212345"
end program elliot8
!end source begin output
299792458.0000000000000000000000000
0.0000012566370614359172885
3.1415926535897931160000000
0.0000000000088541870000000
0.0000000000088541878176204
12345678.0000000000012345678911234
12345678.1234567891123456789212345
I think this version gives as much width on epsilon nought as double
precision can do. I can't believe I got stung with a single-precision
pi. It may have been an over-reaction, but I went through and
appended _dp to every constant. It isn't my idea of pretty, but
wouldn't making everything dp be the best way to avoid having
something implicitly single precision? In particular, it looks odd as
an exponent.
--
Wade Ward
On Jun 4, 9:34 pm, nos@see.signature (Richard Maine) wrote:
> Wade Ward <zaxf@gmail.com> wrote:
> > On May 23, 2:09 am, Jan Vorbrggen <jvorbrueg@not-mediasec.de>
> > wrote:
> > > > it gives the speed of light off by 9 m/s, which is
> > > > faster than Ben Johnson.
> > > Ben Johnson et al. are definitely faster than 9 m/s - they do the 100 m in
> > > less than ten seconds!
> > I realized that right after I posted. I named the only people on the
> > planet who _can_ exceed 9 m/s without machines or free fall or
> > catastrophic consequences.
> Well, it all depends on your frame of reference. Since the speed of
> light was mentioned, seems like relativity ought to be important. I
> figure I'm doing something on the order of 300 m/s sitting here at my
> desk rotating around the center of the earth. Or quite a bit more around
> the Sun or the center of the galaxy. It's enough to make me giddy, or is
> that just from my age? :-)
If we can get creative with frame of reference on running the 100-
meter dash, you could go the Xeno route and simply stay stationary for
the duration of the race. Then define your origin on the sternum of
the fella who crossed the finish line first. You will have gone as
far and as fast as anyone, and the athletes who ran aren't likely to
have a background in special relativity to make the same arguments.
You wouldn't have to worry about a false start. It might work. Then
again, it might not.
--
Wade Ward
|
 |
 |
 |
 |
|