 |
forums.silverfrost.com
Welcome to the Silverfrost forums
|
| View previous topic :: View next topic |
| Author |
Message |
arctica
Joined: 10 Sep 2006
Posts: 105
Location: United Kingdom
|
 Posted: Wed Jul 31, 2013 4:01 pm Post subject: |
 |
|
A test with real data is not working, throwing up a run-time error like before
Floating-point overflow via calls to FOURN and FFTDATA
004027a0 FOURN [+022c] line 264
004022f0 FFTDATA [+00da] line 189 (call fourn)
00401000 main [+0507] line 55 (call fftdata)
The data are a bit over the top in precision so could well be the issue:
76.564239502
76.6045761108
76.5740585327
76.4089889526
76.0467681885
75.442855835
74.5952987671
73.5648574829
72.4396896362
etc
It gets through some of the iteration
iteration 1
rms= 28115.6448491
iteration 2
3.340443934368E+24 9.981365633954E+31
n= 3 Sn/S2 ratio =.547E+04
3.340443934368E+24 3.879393911229E+39
n= 4 Sn/S2 ratio =.341E+08
3.340443934368E+24 1.112922427321E+47
n= 5 Sn/S2 ratio =.183E+12
3.340443934368E+24 3.784955186161E+54
n= 6 Sn/S2 ratio =.106E+16
3.340443934368E+24 6.912957069900E+61
n= 7 Sn/S2 ratio =.455E+19
ress RETURN to close window . . .
at most it should be rounded to 2 decimal places I guess
The main thing is that the program will work, it now seems a case of tweaking the precision settings and/or data.
To view big data files use Notepad++ it does not have the limitation of Excel!
Lester |
|
| Back to top |
|
 |
JohnCampbell
Joined: 16 Feb 2006
Posts: 2554
Location: Sydney
|
| Posted: Thu Aug 01, 2013 1:00 am Post subject: |
|
|
Lester,
Your transform does not look to be converging, by the numbers quoted.
May be this free method doesn't work with real data !!
John |
|
| Back to top |
|
|
arctica
Joined: 10 Sep 2006
Posts: 105
Location: United Kingdom
|
| Posted: Thu Aug 01, 2013 3:12 pm Post subject: |
|
|
Think I will go back and redo the subroutines from scratch rather than try to figure the logic of the published code. With a program that is a computation of a mathematical equation, it is important to document all the key elements. Anyway, it is good practice !
The subroutine FOURN was published in Rick Blakely's book on "Potential Theory and Gravity and Magnetic Applications"
One of the perils of trying to use published code I guess. |
|
| Back to top |
|
|
JohnCampbell
Joined: 16 Feb 2006
Posts: 2554
Location: Sydney
|
| Posted: Fri Aug 02, 2013 1:33 am Post subject: |
|
|
Published ? Did you download it or type it in.
If typed, there are lots of possibilities of typing errors and also mixing characters with bad fonts, eg 1,i,I,l (one, lower case i, upper case I or lower case L) or o/0 or 5/S or 8/B or 9/g. The list could go on.
You could try typing it in twice ( or get someone else to) then use FC or Plato compare to see the result.
John |
|
| Back to top |
|
|
arctica
Joined: 10 Sep 2006
Posts: 105
Location: United Kingdom
|
| Posted: Fri Aug 02, 2013 6:43 am Post subject: |
|
|
It was downloaded code from the journal. I am pretty sure the subroutine FOURN is the same as that shown in the Blakely's book since the variable name declarations are identical. The program code did not cite the original source of the FOURN subroutine code (which they should have done).
Will have another look to see. At least it is a fairly short program!
I do have a Matlab code (again published) which does pretty much the same thing, but looks very cumbersome as the summation is written out in full (n=2 to n=10), so in theory it would not be hard to code it in F95.
I am wondering if there is an issue with the way the existing program is computing the RMS variance for convergence/divergence?
Cheers
Lester |
|
| Back to top |
|
|
JohnCampbell
Joined: 16 Feb 2006
Posts: 2554
Location: Sydney
|
| Posted: Sat Aug 03, 2013 12:14 am Post subject: |
|
|
I think there could be a problem with initialising (zeroing) variables. The following identifies two I think should be done. There are probably many more.
| Code: | subroutine deprms(ht,n,m,rms)
!
real*8 ht(n,m),rms
!
rms = 0 ! initialise rms
!
do i=1,n
do j=1,m
rms=rms+ht(i,j)**2
end do
end do
rms = sqrt ( rms / (n*m) )
return
end
subroutine fftdata (db,data,n,m,n2,n3,nn)
integer*4 n,m,n2,n3
real*8 db(n,m),data(n3)
integer*4 nn(2)
!
integer*4 i,j
!
data = 0 ! initialise data
!
do j=1,m
do i=1,n
data(n2*(j-1)+2*i-1)=db(i,j)
end do
end do
call fourn (data,nn,2,1)
return
end
|
Also I would recommend than NN be declared as integer, rather than real.
Have a look at how I have initialised data and rms. If this should be done, there are probably others like this.
If these changes are required, you would have to question the reliability of the original published code, as a description of rubbish comes to mind.
John |
|
| Back to top |
|
|
arctica
Joined: 10 Sep 2006
Posts: 105
Location: United Kingdom
|
| Posted: Mon Aug 05, 2013 9:36 am Post subject: |
|
|
Thanks for the pointers John. I will look through the code again and check what should and should not be real. There is no ,logic in having the matrix dimension declared real - does seem rather odd coding.
The subroutine FOURN was actually in the Numerical recipes book by Press et al 1986, but don't have a copy to verify that they got the code the same from there.
I managed to get a Matlab code running through Octave to do the same task. However, I will continue to tinker with the Fortran code to see if I can get it functioning properly. Having the thing work ok with test data (synthetics) is all well and good, but the whole idea is to use it "in the field" so to speak.
The code is certainly not well written; more comments are needed for others to follow the logic. Probably a good idea to have an implicit none in to make sure all variables are declared - good programming practice anyway.
Haven't been able to track down the authors yet as the e-mails don't work - always my first port of call.
Will keep you posted.
Cheers
Lester |
|
| Back to top |
|
|
|
|
You cannot post new topics in this forum
You cannot reply to topics in this forum
You cannot edit your posts in this forum
You cannot delete your posts in this forum
You cannot vote in polls in this forum
|
Powered by phpBB © 2001, 2005 phpBB Group
|
FAQ
Search
Memberlist
Usergroups
Register
Profile
Log in to check your private messages
Log in
