Оценить функцию с помощью цикла Fortran 90

Question

Я застрял в процессе, где мне нужно вычислить значения функции f[x,y,z] на сетке. Здесь я рассказываю, как я написал программу, только оценивая на одномерной сетке.

Я написал программу:

program CHISQUARE_MINIMIZATION_VELOCITY_PROFILES 
use distribution 
IMPLICIT none 
integer, parameter :: kp=1001 ! Parameter which states the number of points on the grid. 
integer, parameter :: ndata=13 ! Parameter which states the number of elements of the data file. 
integer, parameter :: nconst=3 ! Fixed integer parameter. 

integer i, j, n 
real*8 rc0, rcf, V00, V0f, d00, d0f, rc, V0, d, z 
real*8 rcr(kp), V0r(kp), d0r(kp), chisq(kp) 

!Scaling radius range 
rc0=0.0d-5  ! kpc 
rcf=1.0d2  ! kpc 
call linspace(rc0,rcf,kp,rcr) 

!**************If I call like this, it works normal***************** 
!CHISQUARED(1.3d0, 130.2d0, 0.12d0, 1.0d0, 1.0d0, 2.0d0, 0.0d0, 0.0d0, 1, !ndata, nconst) 

!  **1.27000000000000  0.745818846396887** 
! Press any key to continue 
!**************If I call like this, it works normal***************** 

    !******* Here is where my problem is**************** 
    do j=1, kp 
    rc=rcr(j) 
    write(*,*) rc, CHISQUARED(rc, 130.2d0, 0.12d0, 1.0d0, 1.0d0, 2.0d0, 0.0d0, 0.0d0, 1, ndata, nconst) 
    enddo 
    !******* Here is where my problem is**************** 

    end program CHISQUARE_MINIMIZATION_VELOCITY_PROFILES 

Я использую модуль, где я вычисление^2 распределения х, исходя из теоретической модели ...

---

MODULE distribution 
IMPLICIT NONE 
CONTAINS 

! I define here the chi^2 function**** 
real*8 function CHISQUARED(rc, V0, d, alpha, gamma, chi, a, b, n, ndata, nconst) 

integer i, n, ndata, nconst 
real*8 rc, V0, d 
real*8 alpha, gamma, chi, a, b, s 
real*8, DIMENSION(ndata,3) :: X 
open(unit=1, file="data.txt") 

s=0.0d0 
do i=1, ndata 
Read(1,*) X(i,:) 
    s=s+((X(i,2)-VELOCITYPROFILE(X(i,1), rc, V0, d, alpha, gamma, chi, a, b, n))/(X(i,3)))**2.0d0 
end do 

CHISQUARED=s/(ndata-nconst) 

end function CHISQUARED 


!****Here I define the model function 
real*8 function VELOCITYPROFILE(r, rc, V0, d, alpha, gamma, chi, a, b, n) 

integer i, n 
real*8 r, rc, V0, d, alpha, gamma, chi, a, b, z 

if (rc < 0.0d0 .OR. d < 0.0d0 .OR. a <0.0d0 .OR. b <0.0d0 .OR. alpha < 0.0d0 .OR. gamma <0.0d0 .OR. chi < 0.0d0 .OR. n<1) then 
    VELOCITYPROFILE=0.0d0 
    return 
    else 
z=0.0d0 
do i=0,n 
    z=z+((V0*((r/rc)**(1.5d0))*(1+a+r/rc)**(-gamma*(2*n+0.5d0)))/((a+(r/rc)**alpha)**(chi/2.0d0)))*(((b+r/rc)**gamma)/d)**i 
end do 

VELOCITYPROFILE=z 
end if 
end function VELOCITYPROFILE 
END MODULE distribution 

** *************** КОНЕЦ МОДУЛЯ ******************************

файл data.txt имеет вид

0.24 37.31 6.15 
0.28 37.92 5.5 
0.46 47.12 3.9 
0.64 53.48 2.8 
0.73 55.14 3.3 
0.82 58.47 2.5 
1.08 66.15 3.3 
1.22 69.39 2.75 
1.45 74.55 5. 
1.71 77.94 2.93 
1.87 81.66 2.5 
2.2  86.81 3.02 
2.28 90.08 2.1 
2.69 94.38 3.92 
2.7  95.36 1.8 

Для того, чтобы получить несколько значений функции CHISQUARED, я использую подпрограмму linspace генерировать разбиение 1-мерной сетки

subroutine linspace(xi,xf,jmax,y) 
integer jmax,j 
real*8 xi,xf,y(jmax) 
y=(/(xi+dble(j-1)*(xf-xi)/(dble(jmax)-1.0d0), j=1, jmax)/) 
end subroutine linspace 

Что происходит, что если в основной программе, Я называю функцию CHISQUARED так:

CHISQUARED(1.3d0, 130.2d0, 0.12d0, 1.0d0, 1.0d0, 2.0d0, 0.0d0, 0.0d0, 1, ndata, nconst) 

    **1.27000000000000  0.745818846396887** 
Press any key to continue 

Я получаю некоторое конечное значение, как, я не знаю, 0,7 или что-то вроде этого. (Я ограничил файл данных, поэтому результат не будет написан, я просто ставлю 0.7 в качестве примера). Однако, когда я положил его в петлю, как это в программе, написанной выше, чтобы получить значения по одной размерной сетке, он дает мне ошибку

**0.000000000000000E+000 NaN** 
forrtl: severe (24): end-of-file during read, unit 1, file C:\Users\Ernesto Lopez Fune\Desktop\Minimize\newone\chisquarerotationcurve\data.txt 
Image    PC  Routine   Line  Source 
chisquarerotation 0040B889 Unknown    Unknown Unknown 
Press any key to continue 

Может кто-нибудь мне посоветовать, что делать в этом случае ? Как преодолеть этот барьер?

Answer 1

Согласно вашей ошибке, вы достигаете конца своего файла. Когда вы вызываете свою подпрограмму один раз, это нормально, но в цикле ваш файл читается несколько раз. После первой итерации ваш файл читается до элемента управления EOF, но для следующей итерации программа больше не может читать, поскольку она уже достигла конца файла.

Вам необходимо использовать инструкцию REWIND(1) до end function CHISQUARED. При этом курсор будет перемещен в начале файла. Кроме того, я думаю, что было бы лучше OPEN ваш файл в основной программе, а не в функции или подпрограмме, чтобы избежать множественного ОТКРЫТИЯ/ЗАКРЫТЬ.

Не забудьте указать CLOSE ваш файл, когда вы закончите иметь дело с ним.