切比雪夫多项式拟合星历
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