song.yz@foxmail.com wechat: math-box

计算数学




切比雪夫多项式拟合星历



ephemeris fit by chebyshev polynomial 

main.app   -----  the fit program  which fit the eph and output the coefficients to file "coe.txt"
main_eval.app --- compute the ephemeris by the chebyshev coefficients which need the "coe.txt"

xyz_ECF ----- input file both in TRS or in TCS

tmp  ----  the temp file for debug 
xyz_com_tmp --- the temp file for debug


src  the source of the program

ephcheby  ----  fit program source
evalcheby  ---- evaluate the coefficients which produced by ephcheby


usage :

       main.app 
       main.app  N tspan
       N is the order of the polynomial 
       tspan is time span for data fit
       
       there also should be a file named xyz_ECF
       
       
       
       main_eval.app
       no input argument
       
       there should be to files 
       coe.txt --- the coefficients file 
       and xyz_ECF  which just read the times of each ephoc and compute the ephemeris

拟合的makefile文件如下。


objects = main.o add_neq.o basecheby.o \
          chol_eq.o chol_rf.o covariance.o inv_dtri.o \
            ephfit.o ephfit_shell.f90

cc = ifort

co = -g -c 
 
main.app:$(objects)
	$(cc) -g -o  ../../main.app $(objects)

main.o: main.f90 
	$(cc) $(co) main.f90

add_neq.o : add_neq.f90
	$(cc) $(co) add_neq.f90

basecheby.o : basecheby.f90 
	$(cc) $(co) basecheby.f90 

chol_eq.o: chol_eq.f90 
	$(cc) $(co) chol_eq.f90
	
chol_rf.o: chol_rf.f90
	$(cc) $(co) chol_rf.f90

covariance.o: covariance.f90
	$(cc) $(co) covariance.f90

inv_dtri.o: inv_dtri.f90
	$(cc) $(co) inv_dtri.f90

ephfit.o: ephfit.f90
	$(cc) $(co) ephfit.f90

ephfit_shell.o: ephfit_shell.f90
	$(cc) $(co) ephfit_shell.f90

subroutine  add_neq(ATA,ATb,oc,H,N)
!---------------------------------subroutine  comment
!  Purpose   :   2018-09-05 13:33           (Created) 
!            add a new observation to normal equations 
!             (ATA+ H'H) x = ATb + H'oc
!  Changes   :
!       *.        
!-----------------------------------------------------
!  Input Parameters    :
!       ATA  ----  initial normal matrix 
!       Atb  ----  initial right vector
!       oc -----   o-c  for liner equation
!       H ------   obs vector
!       N -------  the dimention of the matrix
!  Output Parameters   :    
!      ATA ------  updated normal matrix
!      ATb ------- updated right vector
!-----------------------------------------------------
!  Author       : Song Yezhi        
!  Copyrigt (C) : Shanghai Astronomical Observatory,CAS               
!               (All rights reserved, 2018) 
!----------------------------------------------------
implicit real*8(a-h,o-z)
      
real*8:: ATA(N,N),ATb(N) ,H(N)
      
do i = 1,N
     if (H(i)==0d0) cycle
     do j =1 ,N
       ATA(i,j) = ATA(i,j)+ H(i)*H(j)           
     end do                   
     ATb(i)= ATb(i)+H(i)*oc         
end do 

end subroutine  add_neqsubroutine basecheby(bv,t,N)
!---------------------------------subroutine  comment
!  Version   :  V1.0    
!  Coded by  :  syz 
!  Date      :2010.07.09  
!-----------------------------------------------------
!  Purpose     :  任意阶多项式基函数
! 
!  Post Script :
!       1.
!       2.
!       3.    
!-----------------------------------------------------
!  Input  parameters  :
!       1.    x 标量
!       2.    m 多项式阶数
!  Output parameters  :
!       1.
!       2.
!  Common parameters  :
!       1.
!       2.
!----------------------------------------------------
implicit none

integer        ::N

real*8         ::bv(N+1) ,t 

integer        ::i


bv(1) = 1d0

bv(2) = t 

do i=3,N+1
	
	bv(i) = 2d0* t * bv(i-1)  - bv(i-2)
	
end do 

end subroutine basecheby

subroutine chol_eq(A,b,x,P,N)
!---------------------------------subroutine  comment
!  Version   :  V1.0    
!  Coded by  :  song.yz
!  Date      :  2019.04.15
!-----------------------------------------------------
!  Purpose     : 对称正定方程解算(如法方程)
! 
!  Post Script :
!       1.    修正的不开平方Cholesky分解法
!       2.
!       3.    
!-----------------------------------------------------
!  Input  parameters  :
!       1.  A(N,N)对称正定系数矩阵
!       2.  b(N)右向量
!       3.  N--方程维数
!       4.
!  Output parameters  :
!       1.  x-计算结果
!       2.  P ----协方差 (由于A已经为对称正定矩阵<法矩阵>,因此P为A的逆)
!       3.
!       4.
!-----------------------------------------------------
!  Copyrigt  :   
!        Center for Astro-geodynamics
!        Shanghai Astronomical Observatory
!        Chinese Academy of Sciences   
!----------------------------------------------------
implicit real*8(a-z)

integer::N
real*8::A(N,N),b(n),x(n),P(N,N)

!---------------------------------subroutine variable
real*8::L(N,N),d(n)
integer::i,k

real*8::y(N)

call chol_rf(A,L,d,P,n)
!以上已经得到Cholesky分解

y(1)=b(1)

do i=2,n
   tmp1=0d0
   do k=1,i-1
      tmp1=tmp1+L(i,k)*y(k)
   end do
   y(i)=b(i)-tmp1
end do


x(n)=y(n)/d(n)

do i=n-1,1,-1
  tmp1=0d0
  
  do k=i+1,n
  tmp1=tmp1+L(k,i)*x(k)
  end do
  
  x(i)=y(i)/d(i)-tmp1
end do

end subroutine chol_eq

subroutine chol_rf(A,L,d,P,N)
!---------------------------------subroutine  comment
!  Version   :  V1.0    
!  Coded by  :  syz 
!  Date      :  2019.04.15
!-----------------------------------------------------
!  Purpose     :
! 
!  Post Script :
!       1.  不用开平方的Cholesky分解
!       2.
!       3.    
!-----------------------------------------------------
!  Input  parameters  :
!       1. A(N,N)--输入矩阵
!       2. N----矩阵维书
!       3.
!       4.
!  Output parameters  :
!       1.  L矩阵
!       2.  d---对角矩阵(用向量存储)
!       3.
!       4.  P ---协方差矩阵(由于A已经为对称正定矩阵<法矩阵>,因此P为A的逆)
!-----------------------------------------------------
!  Copyrigt  :   
!        Center for Astro-geodynamics
!        Shanghai Astronomical Observatory
!        Chinese Academy of Sciences   
!----------------------------------------------------
implicit real*8(a-z)
integer::N
real*8::A(n,n),L(n,n),d(n),P(N,N)

!---------------------------------subroutine variable
integer::i,j

real*8::  g(n,n)   

!设置初值
L=0d0

d(1)=a(1,1)

do i=2,n


  do j=1,i-1
   !此层循环算g(i,j) 
     tmp1=0d0
    
     do k=1,j-1
       tmp1=tmp1+g(i,k)*L(j,k)
     end do
     
     g(i,j)=a(i,j)-tmp1
     
  end do
  
  
  do j=1,i-1
  !此层循环算L(i,j)
    L(i,j)=g(i,j)/d(j)
  end do
  
  
  !以下算d(i)
  tmp1=0d0
  do k=1,i-1
   tmp1=tmp1+g(i,k)*L(i,k)
  end do

  d(i)=a(i,i)-tmp1
end do

!设置对角线元素
do i=1,N
   L(i,i)=1d0
end do 


call covariance(L,D,P,N)

end subroutine chol_rf
subroutine covariance(L,D,P,N)
!-----------------------------------------------------
! Purpose : 2019-04-15 13:08                 (Created) 
!  	根据LDL'分解计算协方差矩阵
!              
! Changes :
!            
!-----------------------------------------------------
! Input Parameters    :
!    L   -----
!    D   ---- 对角线元素   
!    N   ---- 矩阵维数   
! Output Parameters   :
!         
!
!-----------------------------------------------------
! Author       : Song Yezhi        
! Copyrigt (C) : Shanghai Astronomical Observatory,CAS               
!               (All rights reserved, 2019) 
!----------------------------------------------------

implicit none

integer     :: N
real*8      :: L(N,N),D(N),P(N,N)

!-----------------------------------------------------
integer     :: i
real*8      :: invL(N,N),invLT(N,N)

!-----------------------------------------------------

 
call inv_dtri(L,invL,N)

invLT = transpose(invL)

do i =1,N
	
	invLT(:,i)=invLT(:,i)/D(i)
	
end do 

P = matmul(invLT,invL)

end subroutine covariance
subroutine ephfit(N,tspan)
!-----------------------------------------------------
! Purpose : 2019-05-08 09:50                 (Created) 
!
!              
! Changes :
!            
!-----------------------------------------------------
! Input Parameters    :
!      M ----- dimention of observation data
!      N ----- order of polynomial  
!       
! Output Parameters   :
!      px ---- coefficients of x    
!      py ---- coefficients of y
!      pz ---- coefficients of z
!-----------------------------------------------------
! Author       : Song Yezhi        
! Copyrigt (C) : Shanghai Astronomical Observatory,CAS               
!               (All rights reserved, 2019) 
!----------------------------------------------------

implicit none


integer             ::i,ioerr,N, M ,iter
 

real*8              :: ATAx(N+1,N+1),ATbx(N+1),px(N+1),covx(N+1,N+1),Hx(N+1)

real*8              :: ATAy(N+1,N+1),ATby(N+1),py(N+1),covy(N+1,N+1),Hy(N+1)

real*8              :: ATAz(N+1,N+1),ATbz(N+1),pz(N+1),covz(N+1,N+1),Hz(N+1)


integer             :: year,mon,day,hour,minu ,day0,hour0,minu0

real*8              :: sec,sec0, x, y ,z  ,t, tspan


open(unit=11,file="xyz_ECF")


open(unit=17,file="tmp")

 

ATAx = 0D0
ATbx = 0d0

ATAy = 0D0
ATby = 0d0

ATAz = 0D0
ATbz = 0d0

iter = 0

do 

   read(11,*,iostat=ioerr) year,mon,day0,hour0,minu0,sec0 , x, y, z
   if(ioerr/=0) exit 
   backspace(11)
   
   
   iter = iter+1
   
   write(12, '(A3,i5)') 'coe',iter
   
   write(17, '(A3,i5)') 'coe',iter
   
   do  
	
	    read(11,*,iostat=ioerr) year,mon,day,hour,minu,sec , x, y, z
	    if(ioerr/=0) exit 
	    
	    
	    
	    !t = day - day0 + hour/24d0 + minu/1440d0 + sec/1440d0/3600d0
	    
	    !t = (day - day0) * 24 - hour0  + hour + minu/60d0 + sec/3600d0
	    
	    t= (day-day0)*24*60 + (hour-hour0)*60 +  (minu - minu0)  + (sec - sec0)/60d0
	    
	    if (t>=tspan) then
	    	backspace(11)
	    	exit
	    end if 
	    
	    t = 2d0*(t-0d0)/tspan - 1d0
	    
	    write(17,'(F15.8,I7,4I4,f10.3,3f20.3)') t, year,mon,day,hour,minu,sec, x, y, z
	    
	    	    
	    x = x/1d3
	    y = y/1d3
	    z = z/1d3
	    
	    call basecheby(Hx,t,N)		
	    call add_NEQ(ATAx,ATbx,x ,Hx,N+1)	
	    
	    call basecheby(Hy,t,N)		
	    call add_NEQ(ATAy,ATby,y ,Hy,N+1)
	    
	    call basecheby(Hz,t,N)		
	    call add_NEQ(ATAz,ATbz,z ,Hz,N+1)
	    
	
   end do 



   call chol_eq(ATAx,ATbx,Px,covx,N+1)

   call chol_eq(ATAy,ATby,Py,covy,N+1)

   call chol_eq(ATAz,ATbz,Pz,covz,N+1)
 
 
   ATAx = 0D0
   ATbx = 0d0

   ATAy = 0D0
   ATby = 0d0

   ATAz = 0D0
   ATbz = 0d0
   
   
   !write(12, '(A3,i5, I7,4I4,f10.3)') 'coe' ,iter , year,mon,day,hour,minu,sec 
   
   write(12,'((F25.12,/))') (px(i),i=1,N+1)
   
   write(12,'((F25.12,/))') (py(i),i=1,N+1)

   write(12,'((F25.12,/))') (pz(i),i=1,N+1)

 

end do 
 

end subroutine ephfit
subroutine ephfit_shell(N,tspan)
!-----------------------------------------------------
! Purpose : 2019-06-05 17:42                 (Created) 
!
!              
! Changes :
!            
!-----------------------------------------------------
! Input Parameters    :
!      
!       
!       
! Output Parameters   :
!         
!
!-----------------------------------------------------
! Author       : Song Yezhi        
! Copyrigt (C) : Shanghai Astronomical Observatory,CAS               
!               (All rights reserved, 2019) 
!----------------------------------------------------

implicit none

integer                  :: N,i

real*8                   :: px(N+1),py(N+1),pz(N+1)

real*8                   :: tspan
!-----------------------------------------------------




open(unit=12,file='coe.txt')


write(12,'(A)') 'order of chebyshev polynomial  and interval of data for fit (minute)'
write(12,'(I5,F18.10)')  N, tspan



call ephfit(N,tspan)


end subroutine ephfit_shell
subroutine inv_dtri(L,invL,N)
!-----------------------------------------------------
! Purpose : 2019-04-15 10:45                 (Created) 
!      计算下三角矩阵的逆矩阵,回带过程,无开平方运算
!              
! Changes :
!            
!-----------------------------------------------------
! Input Parameters    :
!      L ----- 下三角矩阵
!      N ----- 矩阵维数 
!          
! Output Parameters   :
!      invL ------ 下三角的逆矩阵   
!
!-----------------------------------------------------
! Author       : Song Yezhi        
! Copyrigt (C) : Shanghai Astronomical Observatory,CAS               
!               (All rights reserved, 2019) 
!----------------------------------------------------

implicit none

integer   ::   N      
real*8    ::   L(N,N),invL(N,N)

!----------------------------------------------------
integer   ::   i , j ,k 
real*8    ::   xsum
!-----------------------------------------------------


invL = 0d0

invL(1,1) = 1/L(1,1)


do i = 2, N
	
	invL(i,i)= 1/L(i,i)
	
	do j = 1,i-1
		
		xsum = 0d0
		
		do k = 1,i
			xsum = xsum + L(i,k)*invL(k,j)
			
		end do 
		
		invL(i,j)= - xsum / L(i,i)

	end do 
	
end do 


end subroutine inv_dtri
program main
!--------------------------------------------  comment
!  Purpose   :
!  Author    :  Song Yezhi       
!               Song Yezhi  中科院上海天文台
!                       天文地球动力学研究中心
!  Versions and Changes   :
!       V1.0 ----2019-05-08 09:16:37
!
!-----------------------------------------------------
!  Notes     :
!     the coefficients for a polynomial p(x) of degree n that is a best fit (in a least-squares sense) for the data in y. 
!     The coefficients in p are in powers
!
!
!-----------------------------------------------------
!  Copyrigt (C) :
!        Center for Astro-geodynamics
!        Shanghai Astronomical Observatory
!        Chinese Academy of Sciences
!----------------------------------------------------

implicit none

 

integer                            :: N, narg

character*8                        :: strN,strSpan

real*8                             :: tspan    !unit  : minute


narg = iargc()
!get the number of args input 

if (narg<2) then
   write(*,*) 'usage:'
   write(*,*) 'main.app  N  tspan'
   write(*,*) ' N  is the order of the polynomial '
   write(*,*) ' tspan is time span for data fit '
else
  call getarg(1,strN)
  call getarg(2,strSpan)
  
  read(strN,*) N
  read(strSpan,*) tspan 

  call ephfit_shell(N,tspan)
end if 

end program main

用切比雪夫多项式系数,计算星历。 makefile如下



objects = main_eval.o   basecheby.o   chebyval.o  eval.o \
          basevel.o
       
cc = ifort

co = -g -c 
 
main_eval.app:$(objects)
	$(cc) -g -o  ../../main_eval.app $(objects)

 
main_eval.o: main_eval.f90
	$(cc) $(co) main_eval.f90
 
basecheby.o : basecheby.f90 
	$(cc) $(co) basecheby.f90 

 

eval.o: eval.f90
	$(cc) $(co) eval.f90

 
	
chebyval.o: chebyval.f90
	$(cc) $(co) chebyval.f90	
	
basevel.o: basevel.f90
	$(cc) $(co) basevel.f90
 
fortran代码如下。

subroutine basecheby(bv,t,N,bvel)
!---------------------------------subroutine  comment
!  Version   :  V1.0    
!  Coded by  :  syz 
!  Date      :2010.07.09  
!-----------------------------------------------------
!  Purpose     :  切比雪夫多项式基函数
! 
!--------------------------------
!  Input  parameters  :
!       1.    t 自变量
!       2.    m 切比雪夫多项式阶数
!  Output parameters  :
!       1.    bv ----- 切比雪夫基向量
!       2.    bvel ----- 速度基向量(与切比雪夫向量维数相同)
!----------------------------------------------------
implicit none

integer        ::N

real*8         ::bv(N+1) ,t , bvel(N+1)

integer        ::i


bv(1) = 1d0

bv(2) = t 

do i=3,N+1
	
	bv(i) = 2d0* t * bv(i-1)  - bv(i-2)
	
end do 

!位置基向量计算完毕
!---------------------

bvel(1) = 0d0

bvel(2) = 1d0

do i=3, N+1
	
	bvel(i) = 2d0 * bv(i-1) + 2d0 *t * bvel(i-1) - bvel(i-2)
	
end do 

!速度基向量计算完毕


end subroutine basecheby
subroutine basevel(bv,t,N)
!---------------------------------subroutine  comment
!  Version   :  V1.0    
!  Coded by  :  syz 
!  Date      :2010.07.09  
!-----------------------------------------------------
!  Purpose     :  速度基函数
! 
!  Post Script :
!       1.
!       2.
!       3.    
!-----------------------------------------------------
!  Input  parameters  :
!       1.    t 标量
!       2.    N 多项式阶数
!  Output parameters  :
!       1.    bv---速度基向量
!       2.
!----------------------------------------------------
implicit none

integer        ::N

real*8         ::bv(N+1) ,t 

integer        ::i

bv(1) = 0d0         

bv(2) = 1d0
      

do i=3,N 
	
	bv(i) = 2d0*i*t /(i-1d0) * bv(i-1) - i/(i-2d0) * bv(i-2)
	
end do 

end subroutine basevel
subroutine chebyval(p,N,x,fx,vel)
!-----------------------------------------------------
! Purpose : 2019-05-08 10:11                 (Created) 
!    polyval(p,x) returns the value of a polynomial of degree n evaluated at x. 
!    The input argument p is a vector of length n+1 whose elements are the coefficients in
!    powers of the polynomial to be evaluated.           
!    fx = p(1) + p(2)x + p(3)x*x + ....
! Changes :
!            
!-----------------------------------------------------
! Input Parameters    :
!    p  多项式系数
!    N  多项式阶数
!    x  自变量 
! Output Parameters   :
!    fx  应变量    
!
!-----------------------------------------------------
! Author       : Song Yezhi        
! Copyrigt (C) : Shanghai Astronomical Observatory,CAS               
!               (All rights reserved, 2019) 
!----------------------------------------------------

implicit none


integer          :: N ,i
real*8           :: p(N+1),x ,fx ,H(N+1) ,vel, Hvel(N+1)
!-----------------------------------------------------


call basecheby(H,x,N,Hvel)
!获取切比雪夫多项式及其导数的基向量


fx = p(1) * H(1)     ! p(1) *1

vel = p(1) * Hvel(1) ! p(1) *0

do i=2,N+1
	
	fx  = fx + p(i) *H(i)
	
	vel = vel+ p(i)* Hvel(i)

end do 



end subroutine chebyval
subroutine eval(N,tspan)
!-----------------------------------------------------
! Purpose : 2019-05-08 09:50                 (Created) 
!
!              
! Changes :
!            
!-----------------------------------------------------
! Input Parameters    :
!      M ----- dimention of observation data
!      N ----- order of polynomial  
!       
! Output Parameters   :
!      px ---- coefficients of x    
!      py ---- coefficients of y
!      pz ---- coefficients of z
!-----------------------------------------------------
! Author       : Song Yezhi        
! Copyrigt (C) : Shanghai Astronomical Observatory,CAS               
!               (All rights reserved, 2019) 
!----------------------------------------------------

implicit none


integer             :: i,ioerr,N, kkk ,iter ,ibreak
 

real*8              :: px(N+1),py(N+1),pz(N+1)


integer             :: year,mon,day,hour,minu ,day0,hour0,minu0

real*8              :: sec, sec0,fx, fy ,fz  ,t, tspan

character*50        :: line

real*8              :: x,y,z,vx ,vy, vz




open(unit=12,file="xyz_ECF")
          
open(unit=13,file="xyz_com")

open(unit=14,file="xyz_com_tmp")


iter = 0

ibreak = 0

do 
          do 
               read(11,iostat=ioerr,fmt='(A)') line
               if( ioerr/=0) then
               	   ibreak = 1
               	   ! 设置ibreak 主要为跳出外层循环, 也可以直接采用goto 跳出外层循环
               	   exit 
               end if 
               
               ibreak = 0
               
               if (line(1:3)=='coe')  exit

          end do 
          
          
          if (ibreak==1)  exit
          
          
          do kkk = 1,N+1			
          	read(11,*)  px(kkk)
          end do 
          
          do kkk = 1,N+1			
          	read(11,*)  py(kkk)
          end do 	
          
          do kkk = 1,N+1			
          	read(11,*)  pz(kkk)
          end do 	
          
          
          iter = iter+1
          
          
          write(14,'(A3,I5)') 'coe',iter
          
          write(14,'((F25.12,/))') (px(i),i=1,N+1)

          write(14,'((F25.12,/))') (py(i),i=1,N+1)
         
          write(14,'((F25.12,/))') (pz(i),i=1,N+1) 
          
          
          
          read(12,*,iostat=ioerr) year,mon,day0,hour0,minu0,sec0 , x, y, z
          if(ioerr/=0) exit 
          backspace(12)
   
   
           
          do 
          	
          	read(12,*,iostat=ioerr) year,mon,day,hour,minu,sec ,x,y,z   ! ,vx,vy,vz
          	if(ioerr/=0) then
          		ibreak = 1
                exit          		
          	end if 
          	
          	t= (day-day0)*24*60 + (hour-hour0)*60 +  (minu - minu0)  + (sec - sec0)/60d0
	    
	        if (t>=tspan) then
	    	   backspace(12)
	    	   exit
	        end if 
	        
	        t = 2d0*(t-0d0)/tspan - 1d0
          
          	call chebyval(px,N,t,fx,vx)
          	
          	call chebyval(py,N,t,fy,vy)
          	
          	call chebyval(pz,N,t,fz,vz)
          	
          	fx=fx*1d3
          	fy=fy*1d3
          	fz=fz*1d3
          	
          	!----------------------------
          	vx = vx*1d3 /60d0/tspan*2d0
          	vy = vy*1d3 /60d0/tspan*2d0
          	vz = vz*1d3 /60d0/tspan*2d0
          	! it is very important
          	!----------------------------
          	
          	write(13,'(I5,4i4,f10.5,3f20.3,3f20.7,2f15.6,f15.3)')  year,mon,day,hour,minu,sec, fx ,fy,fz,vx ,vy ,vz ,0.0 ,0.0 ,0.0 
          	
          	write(14,'(F15.8,I5,4i4,f10.5,3f20.3,3f20.7,2f15.6,f15.3)') t, year,mon,day,hour,minu,sec, fx,fy,fz,vx ,vy ,vz ,0.0 ,0.0 ,0.0 
          	
          	
          end do 
          
         
          
         
end do 
  

close(12)
close(13)

end subroutine eval
program main_eval
!--------------------------------------------  comment
!  Purpose   :
!  Author    :  Song Yezhi       
!               Song Yezhi  中科院上海天文台
!                       天文地球动力学研究中心
!  Versions and Changes   :
!       V1.0 ----2019-05-08 09:16:37
!
!-----------------------------------------------------
!  Notes     :
!     the coefficients for a polynomial p(x) of degree n that is a best fit (in a least-squares sense) for the data in y. 
!     The coefficients in p are in powers
!
!
!-----------------------------------------------------
!  Copyrigt (C) :
!        Center for Astro-geodynamics
!        Shanghai Astronomical Observatory
!        Chinese Academy of Sciences
!----------------------------------------------------

implicit none

integer           :: N
real*8            :: tspan


open(unit=11,file='coe.txt')

read(11,*)

read(11,*)N ,tspan
 
call eval(N,tspan)


close(11)
end program main_eval