我陷入了需要计算函数值的过程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
我使用计算 chi^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生成一维网格的划分
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
谁能建议我在这种情况下该怎么做？如何克服这个障碍？
