Code debugging

[lehloks]

Hi to all users

I have a data I want to integrate using this function :temp=12.56*(((r(i)+r(1+i)/2+ r(i)/2))**2+((r(i+1)-r(i))*(gr(i)/2+gr(i+1)/2)))
sum_f = sum_f + temp

I have written 2 code to do this,code one is not working properly,I guess its my variables specification. Code 2 its working but only if I integrate a different function which is specified in the code2.Any guidance will help! I have also attached the data I want to integrate.

1

implicit real*4(a-h,o-z)
integer sum_f,i
real r(i),r(i+1)
real,dimension(200) :: r,gr

OPEN(UNIT=10,FILE='RDF.dat',IOSTAT= isost)
i=0
do
i=i+1

read (10,'(F8.6,(5X,F8.6))',end=10) r(i),gr(i)

enddo
10 continue
do i = 1,200
write (*,'(F8.6,(5X,F8.6))') r(i), gr(i)
enddo



r(i)=0.0d0
r(i+1)=1.9750
!n= 5
call trap(r(i),r(i+1),r,gr,sum_f)
write(*,20)r(i),r(i+1),r,gr,sum_f
write(1,20)r(i),r(i+1),r,gr,sum_f
20 format(2x,'r(i)=',f10.5,2x,'r(i+1)=',f10.5,2x,'r=',i4,2x,'gr=',f20.15,'sum_f=',f20.15)

stop
close(1)
end




subroutine trap(r(i),r(i+1),sum_f,r,gr)
!implicit real*4(a-h,o-z)

sum_f = 0

do i,199

temp=12.56*(((r(i)+r(1+i)/2+ r(i)/2))**2+((r(i+1)-r(i))*(gr(i)/2+gr(i+1)/2)))
sum_f = sum_f + temp


end

CODE:2

implicit real*4(a-h,o-z)
integer n,i
real a,b,h,k
real,dimension(200) :: r,gr
OPEN(UNIT=10,FILE='RDF.dat',IOSTAT= isost)
i=0
do
i=i+1
read (10,'(F8.6,(5X,F8.6))',end=10) r(i),gr(i)

enddo
10 continue
do i = 1,200
write (*,'(F8.6,(5X,F8.6))') r(i), gr(i)
enddo



a=0.0d0
b=1.9750
n= 5
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*3.633456E+00*(r**2)*0.041089965


return
end

DATA

1.7750  0.000000E+00
1.8250  0.000000E+00
1.8750  1.088335E-07
1.9250  2.696392E-05
1.9750  7.786679E-04
2.0250  9.117466E-03
2.0750  5.410612E-02
2.1250  1.922786E-01
2.1750  4.879803E-01
2.2250  1.096283E+00
2.2750  2.204320E+00
2.3250  3.438006E+00
2.3750  4.067287E+00
2.4250  3.926939E+00
2.4750  3.366442E+00
2.5250  2.727331E+00
2.5750  2.171983E+00
2.6250  1.735937E+00
2.6750  1.410520E+00
2.7250  1.172828E+00
2.7750  9.996636E-01
2.8250  8.743047E-01
2.8750  7.833833E-01
2.9250  7.172745E-01
2.9750  6.695526E-01
3.0250  6.354343E-01
3.0750  6.120318E-01
3.1250  5.973726E-01
3.1750  5.908303E-01
3.2250  5.906957E-01
3.2750  5.979623E-01
3.3250  6.099009E-01
3.3750  6.279513E-01
3.4250  6.501037E-01
3.4750  6.759086E-01
3.5250  7.047550E-01
3.5750  7.339316E-01
3.6250  7.616673E-01
3.6750  7.856703E-01
3.7250  8.065359E-01
3.7750  8.238194E-01
3.8250  8.399156E-01
3.8750  8.560802E-01
3.9250  8.738620E-01
3.9750  8.954622E-01
4.0250  9.222992E-01
4.0750  9.551873E-01
4.1250  9.967476E-01
4.1750  1.047825E+00
4.2250  1.102791E+00

 

[mikrom]

In your source is lot of errors, for example:

1. as you call subroutine trap you have argument list
r(i), r(i+1), r, gr, sum_f
but in the definition of trap you have argument list
r(i), r(i+1), sum_f, r, gr
The same arguments should be on same place.

2. Why you have the array elements R(i) and R(i+1) in your argument list when you have the whole array R included. IMO the argument list shoud contain the 2 input arrays R and GR and the result sum_f

3. Do you try to sum the array until given length of 199 ? But how do you know that your array has exactly 199 elements ? Array length could be the fourth argument of the subroutine.
This is incorrect too: do i,199.
Please make some effort to study the programming language you want to use...

 

[lehloks]

Hi
I have configured the code a little bit,but I still dont think Im integrating over the data Im reading.Actually I want to have a running integral for the data using the function I have specified or the one I have commented out at the end of the code. Im still working on it.

implicit real*4(a-h,o-z)
integer sum_f,i
!real r(i),r(i+1)
real,dimension(200) :: r,gr


!implicit none
integer, parameter :: n=200
real :: u
!  integer :: i

OPEN(UNIT=10,FILE='RDF.dat',IOSTAT= isost)
i=0
do
i=i+1

read (10,'(F8.6,(5X,F8.6))',end=10) r(i),gr(i)

enddo
10 continue
do i = 1,199
write (*,'(F8.6,(5X,F8.6))') r(i), gr(i)
enddo

call trapezoid_integration(n,5.0)

  contains

    subroutine trapezoid_integration(n,end_val)
      implicit none
      integer :: n
      real ::  end_val
      real :: integral,u,h
      integer :: i

      integral = 0.0

      do i=0,n
         u = (end_val*i)/n

         ! implement Eqn G.4
         if ((i.eq.0).or.(i.eq.n)) then
            integral = integral+integrand(u)
         else
            integral = integral+(2.0*integrand(u))
         end if

      end do

      h=end_val/n
      integral = (h/2.0)*integral

      write (*,*) '#trapezoidal integration = ',integral
    end subroutine trapezoid_integration

    function integrand(x) result (value)
      implicit none
      real :: x
      real :: value

      if (x .lt. 0.00001) then
         x = 0.00001
      end if
!      value = 12.5*(((r(i)+r(1+i)/2+r(i)/2))**2+((r(i+1)-r(i))*(gr(i)/2+gr(i+1)/2))) 
      value = (x**4)*EXP(X)/((EXP(X)-1.0)**2)
    end function integrand
end

 

[lehloks]

Solved.Thanks for the advice!!