Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

[Fortran 90] Output is NaN in variables

  1. Oct 31, 2013 #1
    Hi all,

    For this problem, I have checked all possible factor that may cause NaN in variables such as
    - division by zero --> not found
    - undefined variables - set initialization for variables - ep, hr,ht
    - parameter setting --> double precision problems checked, no issue

    but still couldn't find the cause of the NaN in variable ep, ht and hr as well. need an advice from professionals ( i am an amateur anyway)

    Code (Text):


    !this is code from berenger 2002 for vlf propagation in earth-ionosphere waveguide.
    !equations 7a-7c
    !to date, ionosphere is assumed as free space -29/10/2013

    !**************************************************************************************

    subroutine spherical
    implicit none

    integer                                              :: N, L, i, j, t, k
    integer                                              :: ir !, rr
    double precision                                :: c, f, lamda, dx, dt, sigma, eps0,u0
    double precision                                :: epsr,sigmar,omegap
    double precision                                :: dr, dth,rmh,rph
    double precision,dimension(300,300)   :: ep, ht, hr, ga, gb
    double precision,dimension(300)         :: r
    double precision, parameter               :: pi = 3.14159265d0
    character(len=20)                              :: filename


    !parameter
     k = 201
     L = 100

     f = 3.e9
     u0 = 4*pi*1e-7
     eps0 = 8.8e-12
     c = 1D0/sqrt(u0*eps0)
     lamda = c/f
    ! print *,'c=',c
     
     dx = lamda/20
     dt = dx / (2*c)
     dr = 0.05*lamda
     dth = 1*pi/180          !theta  = pi
     
     !initialize r
     do j = 1,k+1
       r(j) = 0.0
       !print *, 'j=',j,'r(j)=',r(j)
     end do


    !do r = 1,dr*(k-1),dr
     do ir = 1, (k-1)
        r(ir) = (ir - 1)*dr
     end do !IR

     !material parameter , assumed to be free space
      epsr = 1d0
      sigma = 1d0
      sigmar = 1d0
             

    !initialization for updated equations
     do i = 1,k+1
      do j =1 ,k+1
        ep(i,j) = 0
        ht(i,j) = 0
        hr(i,j) = 0
      end do
     end do


    !material properties
    do i = 1, k                                  
        do j = 1,k
            ga(i,j) = exp(-(sigmar*dt)/(eps0*epsr))
            gb(i,j) = (1-exp(-(sigmar*dt)/(eps0*epsr)))/sigmar              
            !print *,'ga(i,j)=',ga(i,j),'gb(i,j)=',gb(i,j)
        end do
    end do

    !**********************************************************************
    !**********************************************************************
     ! update equations
     do t=1,L
     
          write (filename, "('data',I3.3,'.dat')") t
          open (unit=130,file=filename)
               
          ! source
          ep(1,1) = sin(2*pi*f*t*dt)
          !print *, 'ep(1,1)=',ep(1,1)
         
        !  update Ez field
        do i = 2,k
         do j = 2,k
           rph = r(j)+dr/2
           rmh = r(j)-dr/2
           ep(i,j) = ga(i,j)*ep(i,j) + gb(i,j) * (((1/r(j)*dth)*(hr(i,j) - hr(i-1,j)) - (1/r(j)*dr)*(rph*ht(i,j) + rmh*ht(i,j-1))))
           !print*, 'i=',i,'j=',j,'r(j) = ',r(j)
           write (130,*) i,j,ep(i,j)
           if (j == 200) write (130,*) ' '
         end do
        end do
         
         !update Ht (theta) field =Hx
         do  i = 1,k-1
          do j = 1,k-1
             rph = r(j)+dr
             rmh = r(j)
             ht(i,j) = ht(i,j) - (dt/u0/r(j)/dr)*(rph*ep(i,j+1) - rmh*ep(i,j))
             hr(i,j) = hr(i,j) + (dt/u0/r(j)/dth/sin(dth))*(sin(dth)*ep(i+1,j) - sin(dth)*ep(i,j))
             print*,'i=',i,'j=',j,'hr(i,j)  =',hr(i,j)
           end do
          end do
       
       
        close (unit=130)
               
     end do

    end subroutine

     
     
  2. jcsd
  3. Oct 31, 2013 #2

    D H

    User Avatar
    Staff Emeritus
    Science Advisor

    You are dividing by zero (infinity) and then multiplying by zero, which results in NaN.

    Code (Text):

    !initialize r
     do j = 1,k+1
       r(j) = 0.0
       !print *, 'j=',j,'r(j)=',r(j)
     end do

    !do r = 1,dr*(k-1),dr
     do ir = 1, (k-1)
        r(ir) = (ir - 1)*dr
     end do !IR
    The first loop sets the first k+1 elements of r to zero. The second loop sets the first k-1 elements of r to (ir-1)*dr. Elements k and k+1? They're still zero. (Aside: r(1) is also zero.)

    Code (Text):

        !  update Ez field
        do i = 2,k
         do j = 2,k
           rph = r(j)+dr/2
           rmh = r(j)-dr/2
           ep(i,j) = ga(i,j)*ep(i,j) + gb(i,j) * (((1/r(j)*dth)*(hr(i,j) - hr(i-1,j)) - (1/r(j)*dr)*(rph*ht(i,j) + rmh*ht(i,j-1))))
           !print*, 'i=',i,'j=',j,'r(j) = ',r(j)
           write (130,*) i,j,ep(i,j)
           if (j == 200) write (130,*) ' '
         end do
        end do
     
    When j=k that 1/r(j) is 1.0/0.0., or +infinity. On the first pass, hr(i,j) - hr(i-1,j) is 0, and infinity*0 is NaN.
     
  4. Oct 31, 2013 #3
    solved! thank you
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook




Similar Discussions: [Fortran 90] Output is NaN in variables
  1. Fortran 90 (Replies: 1)

  2. Fortran NAN error (Replies: 4)

Loading...