Viewing contents of file '../idllib/contrib/icur/ffcurfit.pro'
;**********************************************************************
pro FFCURFIT,Nfree,weight,FLAMDA,CHISQ1,IFLAG,cii=cii
; called by ffcfit in ffit2
;C**   CURFIT FROM BEVINGTON VIA FWCFIT
;C**   TYPED IN 6/25/82
;C**   modified 7/2/84 to permit freezing of parameters, following
;C**        M. Lampton's WHIZ algorithm
COMMON CFT,X,Y,SIGMAY,eeeeee
COMMON PLCR,YFIT
COMMON CURVE,A,EA,IFIXT
common custompars,dw,lp,x2,x3,x4,x5
common ffits,lu4
if n_elements(lu4) eq 0 then lu4=-1
ea=a*0.
iflag=-1
TW=TOTAL(WEIGHT)
;
n=where(ifixt gt 0,NVAR)
INDX=FIX(N)
;
BETA=dblarr(nvar,/nozero)
ALPHA=dblarr(nvar,nvar,/nozero)
array=dblarr(nvar,nvar,/nozero)
FUN=FUNGUS(x,A,cii=cii)
deriv=GDERIV(x,a)
mm=min(abs(deriv)>0.)
isign=deriv/(abs(deriv)>mm)
logd=alog10(abs(deriv))
logw=alog10(weight+1.E-9)
z=weight*(y-fun)
for j=0,nvar-1 do beta(j)=total(z*deriv(*,indx(j)))
;
for j=0,nvar-1 do begin            
   JI=INDX(J)
   for k=0,j DO begin 
      KI=INDX(K)
      z=(logW+logD(*,JI)+logD(*,KI))>(-37.)
      ALPHA(J,K)=TOTAL(isign(*,ji)*isign(*,KI)*10^(z)) 
      endfor
   endfor
;
yfit=fungus(x,a,cii=cii)
IFLAG=-1
for J=0,NVAR-1 do begin
   for K=0,J do ALPHA(K,J)=ALPHA(J,K)
   IF ALPHA(J,J) EQ 0. then goto,curerr
   endfor
IFLAG=0
CHISQ1=FCHISQ(NFREE,y,yfit,weight)
lab71:
sqralpha=sqrt(ABS(alpha))
for J=0,NVAR-1 do begin
   dj=sqrt(alpha(j,j))
   for K=0,NVAR-1 do ARRAY(J,K)=ALPHA(J,K)/dj/sqrALPHA(K,K)
   ARRAY(J,J)=1.D0+FLAMDA
   endfor
;
iflag=-99
arr=array
array=invert(arr)
;
b=a
dal=dblarr(nvar,/nozero)
for j=0,nvar-1 do dal(j)=alpha(j,j)
dal=sqrt(dal)
for J=0,NVAR-1 do begin
   JI=INDX(J)
   b(ji)=b(ji)+total(beta*array(j,*)/dal(j)/dal)
;for K=0,NVAR-1 do B(JI)=b(ji)+BETA(K)*ARRAY(J,K)/SQRT(ALPHA(J,J))/sqrt(ALPHA(K,K))
   endfor
YFIT=FUNGUS(x,B,cii=cii)
CHISQR=FCHISQ(NFREE,y,yfit,weight)
case 1 of
   chisq1 lt chisqr: begin
      FLAMDA=10.D0*FLAMDA
      IF FLAMDA LT 101.D0 then GOTO,lab71
      IFLAG=-92
      RETURN
      end
   chisq1 ge chisqr: begin
      a=b
      saa=SQRT(ABS(ARRAY/ALPHA))
      for J=0,NVAR-1 do begin
         JI=INDX(J)
;         EA(JI)=SQRT(ARRAY(J,J)/ALPHA(J,J))
         EA(JI)=SAA(J,J)
         endfor
      end
   endcase
FLAMDA=FLAMDA/10.D0
CHISQ1=CHISQR
iflag=0
RETURN
;
curerr: 
bell
print,' *** CURFIT failed to converge'
if lu4 le 0 then return
printf,lu4,'*  CURFIT ERRORS: IFLAG,J,Nvar=',IFLAG,J,Nvar
printf,lu4,'*  I,A(I),DERIV(I),BETA(I),ALPHA(I,J=(1,N))'
for I=0,Nvar-1 do begin
   ji=indx(i)
   printf,lu4,'* ',I,A(jI),DERIV(jI),BETA(I),(ALPHA(I,0:nvar-1))
   endfor
RETURN
END