Intergrating data in an external folder 1

[lehloks]

Hi everyone
I am tryying to integrate data in an external file using trapezium rule. I can compile the code and run it. But it does not seem to intergrate data in an external folder although i get the output. It looks like its intergrating the function and its parameters only. any help will be much appreciated.
I have attached the code and the data I want to intergrate below:

regards
Lehloks
CODE:
implicit real*4(a-h,o-z)
integer n,ok
real a,b,h,gr,r
logical ex
!real,allocatable,dimension) ::gr, r

! ** OPEN DATA FILE AND RESULTS FILE **

!OPEN(UNIT=10,FILE='RDF.dat',STATUS='OLD')
!OPEN(UNIT=10,FILE='RDF.dat',IOSTAT= ok)

!read(10,'(i1)')n
!allocate(gr)
!allocate(r)
!inquire(file='RDF.dat', exist=ex)
!inquire(file=name, exist=ex)
!if (ex) then
OPEN(UNIT=10,FILE='RDF.dat',IOSTAT= isost)
do
!inquire(file='RDF.dat', exist=ex)
read (10,'(F8.6,(5X,F8.6))',end=10) gr,r
write(*,*) iost
!write (*,'(F813.6,(5X,F8.6))') gr, r
enddo

10 continue



a=0.0d0
b=1.0d0
n=90
ok=0
call trap(a,b,n,apot)
write(*,20)a,b,n,apot
write(1,20)a,b,n,apot
20 format(2x,'a=',f10.5,2x,'b=',f10.5,2x,'n=',i4,2x,'apot=',f20.15)

stop
close(1)
end




subroutine trap(a,b,n,s)
implicit real*4(a-h,o-z)
f(r) = 4*3.1415927*0.2*(r**2)*2
h=(b-a)/dfloat
s=f(a)
r=a
nm=n/2
do 30 i=1,nm
x=x+h
30 s=s/2+f(r)
s=s*h
return
end

Data:


0.025000 0.000000
0.075000 0.000000
0.125000 0.000000
0.175000 0.000000
0.225000 0.000000
0.275000 0.000000
0.325000 0.000000
0.375000 0.000000
0.425000 0.000000
0.475000 0.000000
0.525000 0.000000
0.575000 0.000000
0.625000 0.000000
0.675000 0.000000
0.725000 0.000000
0.775000 0.000000
0.825000 0.000000
0.875000 0.000000
0.925000 0.000000
0.975000 0.000000
1.025000 0.000000
1.075000 0.000000
1.125000 0.000000
1.175000 0.000000
1.225000 0.000000
1.275000 0.000000
1.325000 0.000000
1.375000 0.000000
1.425000 0.000000
1.475000 0.000000
1.525000 0.000000
1.575000 0.000000
1.625000 0.000000
1.675000 0.000000
1.725000 0.000000
1.775000 0.000000
1.825000 0.000000
1.875000 0.000020
1.925000 0.000091
1.975000 0.000354
2.025000 0.001007
2.075000 0.002950
2.125000 0.007610
2.175000 0.018784
2.225000 0.041836
2.275000 0.089604
2.325000 0.231973
2.375000 0.673571
2.425000 1.460788
2.475000 2.037844
2.525000 1.961896
2.575000 1.686782
2.625000 1.777382
2.675000 2.267358
2.725000 2.761957
2.775000 2.934413
2.825000 2.957162
2.875000 3.194042
2.925000 3.791247
2.975000 4.695648
3.025000 5.840049
3.075000 7.093918
3.125000 8.050863
3.175000 8.157009
3.225000 7.541305
3.275000 6.664297
3.325000 5.612486
3.375000 4.408674
3.425000 3.404541
3.475000 2.945315
3.525000 2.637840
3.575000 2.121458
3.625000 1.562140
3.675000 1.225294
3.725000 1.067429
3.775000 1.007232
3.825000 1.010981
3.875000 1.084994
3.925000 1.209127
3.975000 1.358738
4.025000 1.522696
4.075000 1.673802
4.125000 1.818884
4.175000 1.968916
4.225000 2.123943
4.275000 2.285839
4.325000 2.439340
4.375000 2.530822
4.425000 2.514522
4.475000 2.371838
4.525000 2.157134
4.575000 1.948307
4.625000 1.782963
4.675000 1.689690
4.725000 1.667512
4.775000 1.704558
4.825000 1.786024
4.875000 1.888476
4.925000 1.991947
4.975000 2.074554
5.025000 2.125446
5.075000 2.143362
5.125000 2.156244
5.175000 2.201955
5.225000 2.311914
5.275000 2.495983
5.325000 2.726427
5.375000 2.904423
5.425000 2.911576
5.475000 2.747076
5.525000 2.538788
5.575000 2.409045
5.625000 2.418356
5.675000 2.564244
5.725000 2.783672
5.775000 3.017873
5.825000 3.198384
5.875000 3.273338
5.925000 3.176226
5.975000 2.908067
6.025000 2.576140
6.075000 2.280277
6.125000 2.069323
6.175000 1.913700
6.225000 1.804934
6.275000 1.738271
6.325000 1.712798
6.375000 1.720140
6.425000 1.742850
6.475000 1.773337
6.525000 1.808804
6.575000 1.861121
6.625000 1.918463
6.675000 1.973188
6.725000 2.008807
6.775000 2.004753
6.825000 1.955524
6.875000 1.898719
6.925000 1.875917
6.975000 1.918234
7.025000 2.013022
7.075000 2.116971
7.125000 2.200793
7.175000 2.257111
7.225000 2.317633
7.275000 2.392609
7.325000 2.476039
7.375000 2.551951
7.425000 2.599623
7.475000 2.608005
7.525000 2.579772
7.575000 2.524388
7.625000 2.470380
7.675000 2.406537
7.725000 2.323634
7.775000 2.214267
7.825000 2.090324
7.875000 1.975732
7.925000 1.886941
7.975000 1.842203
8.025000 1.837723
8.075000 1.876098
8.125000 1.945028
8.175000 2.019817
8.225000 2.068400
8.275000 2.094981
8.325000 2.100951
8.375000 2.103055
8.425000 2.088477
8.475000 2.045760
8.525000 1.972763
8.575000 1.884686
8.625000 1.808937
8.675000 1.754546
8.725000 1.726358
8.775000 1.724135
8.825000 1.738965
8.875000 1.756061
8.925000 1.753199
8.975000 1.717060
9.025000 1.658192
9.075000 1.601135
9.125000 1.560162
9.175000 1.556259
9.225000 1.587308
9.275000 1.638478
9.325000 1.688446
9.375000 1.720518
9.425000 1.732540
9.475000 1.725602
9.525000 1.704314
9.575000 1.687294
9.625000 1.677712
9.675000 1.677202
9.725000 1.683873
9.775000 1.687297
9.825000 1.677633
9.875000 1.663318
9.925000 1.647225
9.975000 1.639291

 

[FJacq]

Sorry but I don't understand your programming : you read the data file but you don't use what you read (the variable gr and r are never used after having been set by the READ statement).

François Jacq

 

[lehloks]

Dear François Jacq

Im not familier with fortran but I have included my gr and r in my function at the end.If there is anyway can you please guid me on how can I go about doing this.

Regards
Lehloks

 

[lehloks]

Dear François Jacq subroutine trap(a,b,n,s)
I have changed x = x + h into r = r + h
gr is suppose to be a constant in the code.


implicit real*4(a-h,o-z)
f(r) = 4*3.1415927*0.2*(r**2)*2
h=(b-a)/dfloat
s=f(a)
r=a
nm=n/2
do 30 i=1,nm
r=r+h

30 s=s/2+f(r)
s=s*h
return
end

 

[salgerman]

Actually, I don't even see where the entire file is being read...the do loop reads scalars gr and r over and over and over, which means that at the end, only the values from the last row persist.

You need declare those variable arrays and read a row per index value inside the loop.

 

[FJacq]

You don't understand how a programming language works. This is not related to Fortran only because C, C++ JAVA, Pascal ... work in the same way.

Look into your function subroutine TRAP : it receives, as arguments, the variables a, b , n, s and that's all. This subroutine also uses the symbols f, r, h, s, nm ... which are all INTERNAL definitions. In particular, the variable r of TRAP has nothing to do with the variable r of your main program.

François Jacq

 

[mikrom]

I don't understand right what your code posted above should mean, but:

IMO, you have given a tabulated function y = f(x):

x[sub]0[/sub], y[sub]0[/sub] 
...
x[sub]k-1[/sub], y[sub]k-1[/sub] 
x[sub]k[/sub], y[sub]k[/sub] 
...
x[sub]N[/sub], y[sub]N[/sub]

and you need to integrate it, i.e. to compute numerical quadrature using tarpezoid rule.

If so - it looks to be simple: you have to sum surface areas of all elementary trapeziums.
If I'm not wrong, it must be something like this:

sum = 0
foreach k in 1,..,N do:
  sum = sum + (x[sub]k+1[/sub] - x[sub]k[/sub])*(y[sub]k[/sub] + y[sub]k+1[/sub])/2

 

[mikrom]

but truly said, it looks for me like homework...

 

[mikrom]

Here is what I proposed:

[COLOR=#0000ff]! trapezoidal rule for integrating dataset[/color]
[COLOR=#0000ff]! x1, y1[/color]
[COLOR=#0000ff]! ...[/color]
[COLOR=#0000ff]! x_n, y_n[/color]
[COLOR=#a020f0]program[/color] integrate_dataset
  [COLOR=#2e8b57][b]implicit[/b][/color] [COLOR=#2e8b57][b]none[/b][/color]
[COLOR=#2e8b57][b]  real[/b][/color]:: x, y, x_old, y_old, sum_trapez
  [COLOR=#2e8b57][b]integer[/b][/color] :: num_line, stat

  num_line [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]0[/color]
  sum_trapez [COLOR=#804040][b]=[/b][/color] [COLOR=#ff00ff]0[/color]
  [COLOR=#804040][b]do[/b][/color]
    [COLOR=#804040][b]read[/b][/color]([COLOR=#804040][b]*[/b][/color], [COLOR=#804040][b]*[/b][/color], [COLOR=#804040][b]end[/b][/color][COLOR=#804040][b]=[/b][/color][COLOR=#804040][b]90[/b][/color], [COLOR=#804040][b]iostat[/b][/color][COLOR=#804040][b]=[/b][/color]stat) x, y
    [COLOR=#804040][b]if[/b][/color] (stat [COLOR=#804040][b].ne.[/b][/color] [COLOR=#ff00ff]0[/color]) [COLOR=#804040][b]then[/b][/color]
      [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'* Error reading data on line: '[/color], num_line
      [COLOR=#804040][b]go to[/b][/color] [COLOR=#ff00ff]99[/color]
    [COLOR=#804040][b]end if[/b][/color]
    num_line [COLOR=#804040][b]=[/b][/color] num_line [COLOR=#804040][b]+[/b][/color] [COLOR=#ff00ff]1[/color]
    [COLOR=#804040][b]if[/b][/color] (num_line [COLOR=#804040][b]>[/b][/color] [COLOR=#ff00ff]1[/color]) [COLOR=#804040][b]then[/b][/color]
      sum_trapez [COLOR=#804040][b]=[/b][/color] sum_trapez [COLOR=#804040][b]+[/b][/color] (x [COLOR=#804040][b]-[/b][/color] x_old)[COLOR=#804040][b]*[/b][/color](y [COLOR=#804040][b]+[/b][/color] y_old)[COLOR=#804040][b]/[/b][/color][COLOR=#ff00ff]2[/color]
    [COLOR=#804040][b]end if[/b][/color]
    [COLOR=#0000ff]!write(*, *) 'current line: ', x, y, 'previous ', x_old, y_old, &[/color]
    [COLOR=#0000ff]!            'sum: ', sum_trapez[/color]
    [COLOR=#0000ff]!    [/color]
    x_old [COLOR=#804040][b]=[/b][/color] x
    y_old [COLOR=#804040][b]=[/b][/color] y
  [COLOR=#804040][b]end do[/b][/color]

[COLOR=#804040][b]  90 continue[/b][/color]
  [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'Number of data lines processed: '[/color], num_line
  [COLOR=#804040][b]write[/b][/color]([COLOR=#804040][b]*[/b][/color],[COLOR=#804040][b]*[/b][/color]) [COLOR=#ff00ff]'Computed quadrature using trapezouid rule S = '[/color], sum_trapez  
  [COLOR=#0000ff]! [/color]
[COLOR=#804040][b]  99 continue[/b][/color]
[COLOR=#a020f0]end program[/color] integrate_dataset

Output:

$ g95 integrate_dataset.f95 -o integrate_dataset

$ integrate_dataset < integrate_dataset.dat
 Number of data lines processed:  200
 Computed quadrature using trapezouid rule S =  17.718332

 

[lehloks]

mikrom

Thanks a lot you are a star. Its working properly.