perm filename CURFIT.F4[NET,GUE] blob sn#041474 filedate 1973-05-14 generic text, type T, neo UTF8
      subroutine curfit(x,y,sigmay,npts,nterms,mode,a,deltaa,
     c sigmaa,flamda,yfit,chisqr)
      implicit double precision (a-h,o-z)
      dimension x(100),y(100),sigmay(100),a(10),deltaa(10),sigmaa(10),
     c yfit(100)
      dimension weight(100),alpha(10,10),beta(10),deriv(10),
     c array(10,10),b(10)
11    nfree=npts-nterms
      if(nfree) 13,13,20
13    chisqr=0.
      goto 110
20    do 30 i=1,npts
21    if (mode) 22,27,29
22    if(y(i)) 25,27,23
23    weight(i)=1./y(i)
      goto 30
25    weight(i)=1./(-y(i))
      goto 30
27    weight(i)=1.
      goto 30
29    weight(i)=1./sigmay(i)**2
30    continue
31    do 34 j=1,nterms
      beta(j)=0.
      do 34 k=1,j
34    alpha(j,k)=0.
41    do 50 i=1,npts
      call fderiv(x,i,a,deltaa,nterms,deriv)
      do 46 j=1,nterms
      beta(j)=beta(j)+weight(i)*(y(i)-functn(x,i,a))*deriv(j)
      do 46 k=1,j
46    alpha(j,k)=alpha(j,k)+weight(i)*deriv(j)*deriv(k)
50    continue
51    do 53 j=1,nterms
      do 53 k=1,j
53    alpha(k,j)=alpha(j,k)
61    do 62 i=1,npts
62    yfit(i)=functn(x,i,a)
63    chisq1=fchisq(y,sigmay,npts,nfree,mode,yfit)
71    do 74 j=1,nterms
      do 73 k=1,nterms
73    array(j,k)=alpha(j,k)/dsqrt(alpha(j,j)*alpha(k,k))
74    array(j,j)=1.+flamda
80    call matinv(array,nterms,det)
81    do 84 j=1,nterms
      b(j)=a(j)
      do 84 k=1,nterms
84    b(j)=b(j)+beta(k)*array(j,k)/dsqrt(alpha(j,j)*alpha(k,k))
91    do 92 i=1,npts
92    yfit(i)=functn(x,i,b)
93    chisqr=fchisq(y,sigmay,npts,nfree,mode,yfit)
      if (chisq1-chisqr) 95,101,101
95    flamda=10.*flamda
      goto 71
101   do 103 j=1,nterms
      a(j)=b(j)
103   sigmaa(j)=dsqrt(array(j,j)/alpha(j,j))
      flamda=flamda/10.
110   return
      end