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

计算数学



数据拟合



多项式拟合星历

本节与下一节内容介绍用多项式和切比雪夫多项式拟合星历问题。实际计算中,由于数据较大,一般通过 逐行读取文件的方法。本节代码可以编译成两个程序,一个是拟合程序,生成拟合系数。一个是通过系数计算星历程序。 下一节内容用法同本节。
拟合的Makefile文件如下。

objects = main.o add_neq.o basefunc.o \
          chol_eq.o chol_rf.o covariance.o inv_dtri.o \
          polyfit.o polyval.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

basefunc.o : basefunc.f90 
	$(cc) $(co) basefunc.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

polyfit.o: polyfit.f90
	$(cc) $(co) polyfit.f90
	
polyval.o: polyval.f90
	$(cc) $(co) polyval.f90	
	
ephfit_shell.o: ephfit_shell.f90
	$(cc) $(co) ephfit_shell.f90
	
利用make规则,还可以写成。

objects = main.obj add_neq.obj basefunc.obj \
          chol_eq.obj chol_rf.obj covariance.obj inv_dtri.obj \
          polyfit.obj polyval.obj  ephfit.obj ephfit_shell.obj

cc = ifort

co = -c 

main.exe:$(objects)
	$(cc) -o  main.exe $(objects)

$(obj): %.obj: %.f90
	$(cc) $(co) $< $@
拟合算法fortran代码如下。

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_neq
subroutine basefunc(bv,x,m)
!---------------------------------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 real*8(a-z)

integer::M

real*8::bv(m)

integer::i

bv(1)=1d0

do  i=2,m
    bv(i)=bv(i-1)*x
end do
end subroutine basefunc
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 
	    
	    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 basefunc(Hx,t,N+1)		
	    call add_NEQ(ATAx,ATbx,x ,Hx,N+1)	
	    
	    call basefunc(Hy,t,N+1)		
	    call add_NEQ(ATAy,ATby,y ,Hy,N+1)
	    
	    call basefunc(Hz,t,N+1)		
	    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 polynomial  and end of time'
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       
!              
!  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
subroutine  polyfit(x,y,M,p,N)
!-----------------------------------------------------
! Purpose : 2019-05-08 09:40                 (Created) 
!    
!    The coefficients in p (px =p(1) + p(2)x + p(3)x*x + .... )          
! Changes :
!            
!-----------------------------------------------------
! Input Parameters    :
!      x  ----  自变量数组
!      y  ----  应变量数组
!      M  ----- 数据维数
!      N -----  多项式阶数 
! Output Parameters   :
!     p --参数    P维数为多项式阶数+1
!
!-----------------------------------------------------
! Author       : Song Yezhi        
! Copyrigt (C) : Shanghai Astronomical Observatory,CAS               
!               (All rights reserved, 2019) 
!----------------------------------------------------

implicit none

integer          :: N,M
real*8           :: x(M),y(M),p(N+1)
integer          :: i
real*8           :: H(N+1),ATA(N+1,N+1),ATb(N+1),PX(N+1,N+1)

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

ATA =0D0
ATb =0d0

do i=1,M
	call basefunc(H,x(i),N+1)
	
	
	call add_NEQ(ATA,ATb,y(i),H,N+1)	
	
end do 
    
call chol_eq(ATA,ATb,Px,p,N+1)


end subroutine  polyfitsubroutine polyval(p,N,x,fx)
!-----------------------------------------------------
! 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)
!-----------------------------------------------------
call basefunc(H,x,N+1)
fx =p(1)
do i=2,N+1	
	fx =fx + p(i) *H(i)
end do 
end subroutine polyval
利用系数计算星历的makefile如下。

objects = main_eval.o   basefunc.o   polyval.o  eval.o \
          polyvel.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
 
basefunc.o : basefunc.f90 
	$(cc) $(co) basefunc.f90 

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

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

polyvel.o: polyvel.f90
	$(cc) $(co) polyvel.f90
利用系数计算星历的fortran代码如下。

subroutine basefunc(bv,x,m)
!---------------------------------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 real*8(a-z)

integer::M

real*8::bv(m)

integer::i

bv(1)=1d0

do  i=2,m
    bv(i)=bv(i-1)*x
end do
end subroutine basefunc
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  
          	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 
          
          	call polyval(px,N,t,fx)
          	
          	call polyval(py,N,t,fy)
          	
          	call polyval(pz,N,t,fz)
          	
          	call polyvel(px,N,t,vx)
          	
          	call polyvel(px,N,t,vy)
          	
          	call polyvel(pz,N,t,vz)
          	
          	fx=fx*1d3
          	fy=fy*1d3
          	fz=fz*1d3
          	
          	!-------------------------
          	vx = vx*1d3 /60d0 
          	vy = vy*1d3 /60d0 
          	vz = vz*1d3 /60d0 
          	!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
subroutine polyval(p,N,x,fx)
!-----------------------------------------------------
! 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)
!-----------------------------------------------------


call basefunc(H,x,N+1)


fx =p(1)

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

end do 

end subroutine polyval
subroutine polyvel(p,N,t,vel)
!-----------------------------------------------------
! Purpose : 2019-06-05 19:24                 (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              :: p(N+1)

real*8              :: vel ,t


real*8              :: H(N+1)

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

H(1) = 0d0

H(2) = 1d0

H(3) = t

do i = 4,N+1
	
    H(i) = H(i-1) *t
	
end do 
!这里H仅仅包含多项式,不包含前面导数系数,导数系数放在下面循环计算
vel = 0d0

do i= 1,N+1
	
	vel = vel + H(i)*p(i) * (i-1)
	
end do 
 
end subroutine polyvel