40 lines
1.1 KiB
FortranFixed
40 lines
1.1 KiB
FortranFixed
SUBROUTINE svdfit(x,y,sig,ndata,a,ma,u,v,w,mp,np,chisq,funcs)
|
|
INTEGER ma,mp,ndata,np,NMAX,MMAX
|
|
double precision chisq,a(ma),sig(ndata),u(mp,np),v(np,np),w(np),
|
|
*x(ndata),y(ndata),TOL
|
|
c mp>=ndata, np>=ma. ma is the number of coefficients
|
|
EXTERNAL funcs
|
|
PARAMETER (NMAX=1000,MMAX=50,TOL=1.0d-10)
|
|
CU USES svbksb,svdcmp
|
|
INTEGER i,j
|
|
double precision sumup,thresh,tmp,wmax,afunc(MMAX),b(NMAX)
|
|
do 12 i=1,ndata
|
|
call funcs(x(i),afunc,ma,i)
|
|
tmp=1.0d0/sig(i)
|
|
do 11 j=1,ma
|
|
u(i,j)=afunc(j)*tmp
|
|
11 continue
|
|
b(i)=y(i)*tmp
|
|
12 continue
|
|
call svdcmp(u,ndata,ma,mp,np,w,v)
|
|
wmax=0.0d0
|
|
do 13 j=1,ma
|
|
if(w(j).gt.wmax)wmax=w(j)
|
|
13 continue
|
|
thresh=TOL*wmax
|
|
do 14 j=1,ma
|
|
if(w(j).lt.thresh)w(j)=0.0d0
|
|
14 continue
|
|
call svbksb(u,w,v,ndata,ma,mp,np,b,a)
|
|
chisq=0.0d0
|
|
do 16 i=1,ndata
|
|
call funcs(x(i),afunc,ma,i)
|
|
sumup=0.0d0
|
|
do 15 j=1,ma
|
|
sumup=sumup+a(j)*afunc(j)
|
|
15 continue
|
|
chisq=chisq+((y(i)-sumup)/sig(i))**2
|
|
16 continue
|
|
return
|
|
END
|