Numerical errors in differentiation program

[BerkeleyEric]

Hi all,

I am new to Fortan and am having some trouble with the first program I wrote. I have a piece of software that outputs a large file of formatted data (heat.dat) that looks like:

# Comment line
1 -8082.771321 -8080.455672 -7850.944265
2 -8082.861127 -8080.528031 -7851.010982
3 -8082.870879 -8080.515755 -7851.007309
4 -8082.880682 -8080.50318 -7851.003825
...

First there is a comment line. The subsequent lines contain four columns: a "step" number and three values. Each step corresponds to a fixed amount of time.

I need to take the time-derivative of the values, and then multiply by a fixed constant to get things in the right units.

Here is the program I wrote:

IMPLICIT NONE
INTEGER :: OpenStatus, Counter, TotalSteps, Discard
REAL :: dRx, dRy, dRz, Rx1, Ry1, Rz1, Rx2, Ry2, Rz2, TimeStep_ps, TimeStep_s, Convert


! Input the heat.dat file
OPEN (UNIT = 1, FILE = "heat.dat", STATUS = "OLD", ACTION = "READ", POSITION = "REWIND", IOSTAT = OpenStatus)

! Give an error if the file does not open correctly
IF (OpenStatus > 0) STOP "*** Cannot open heat.dat ***"

! Prompt the user for the simulation time step
PRINT *, "What is the time step in units of picoseconds?"
READ *, TimeStep_ps
TimeStep_s = TimeStep_ps * 1E-12

! Determine the total number of time steps in the heat.dat file and move to the first line of data
Counter = 0
DO
READ(1, *, end=100, err=100)
Counter = Counter + 1
END DO
100 CONTINUE
TotalSteps = Counter
REWIND 1
READ (1, *)


! Create the deriv.dat file
OPEN (UNIT = 2, FILE = "deriv.dat", STATUS = "REPLACE", ACTION = "WRITE", IOSTAT = OpenStatus)

! Give an error if the file does not open correctly
IF (OpenStatus > 0) STOP "*** Cannot open deriv.dat ***"

! Set the proper conversion factor to get units of erg*cm
Convert = 1.60217646D-20

! Loop to calculate the derivatives and write the output
DO Counter = 1, TotalSteps - 1, 1

READ (1, *) Discard, Rx1, Ry1, Rz1
READ (1, *) Discard, Rx2, Ry2, Rz2

dRx = ((Rx2 - Rx1)/TimeStep_s)*Convert
dRy = ((Ry2 - Ry1)/TimeStep_s)*Convert
dRz = ((Rz2 - Rz1)/TimeStep_s)*Convert

WRITE (2, fmt=800) Counter, dRx, dRy, dRz
800 FORMAT (I9, 1X, E30.20, 1X, E30.20, 1X, E30.20)

! Back up to the beginning of the preceding line in heat.dat
BACKSPACE 1

END DO

! Close the files heat.dat and deriv.dat
CLOSE (1)
CLOSE (2)


END PROGRAM HeatDeriv


It seems like the program is working, but the numbers don't seem to be exactly right. After a couple of decimal points the numbers seem wrong (at least compared to trying some in my graphing calculator). I read a little on numerical errors from conversion between normal numbers and binary, but I am not positive this is the problem and am unsure how to fix this.

Thanks for any help you can provide!
Eric

 

[mikrom]

Look at this thread how to compute derivatives from a data set:

http://www.tek-tips.com/viewthread.cfm?qid=1567028

 

[BerkeleyEric]

mikrom,

Thanks. I think I fixed the problem by declaring the variables with "DOUBLE PRECISION".

Cheers,
Eric