Initial commit

1 files changed, 400 insertions, 0 deletions

diff --git a/noao/onedspec/fortran/nlcfit.f b/noao/onedspec/fortran/nlcfit.f
new file mode 100644
index 00000000..80aff616
--- /dev/null
+++ b/noao/onedspec/fortran/nlcfit.f
@@ -0,0 +1,400 @@
+      SUBROUTINE NLCFIT(IM,INN,IN,INTA,XEPSI,XV,XYD)
+C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+C     NONLINEAR LEAST SQUARES FITTING USING SIMPLEX
+C     METHOD AND QUADRATIC APPROXIMATION.
+C     WITH LINEAR PARAMETER ELIMINATION.
+C-------------------------------------------------------------
+      INTEGER IM,INN,IN,INTA
+      REAL XEPSI,XV(IM),XYD(IM)
+      COMMON /NLC/ EPSI,IFLAG,IL,IQ,INDEX,F(15,120),M,N
+      COMMON /NLC/ SOLD,Y(20),YVAL,XF(11),X(11,11),V(120),YD(120,10)
+      COMMON /NLC/ GG(11,11),GINV(11,11),EM(15,120),BB(5),NT
+      COMMON /NLC/ GA(120,5),NN
+	COMMON /NLCOUT/ FF(120),PARS(10),BPARS(10)
+      DIMENSION SUMC(11),XC(11),XE(11),XCO(11),XR(11)
+      REAL LERROR
+C-----
+C RESET ERROR HANDLER
+c...UNIX has general handler only!
+c	call trpfpe (0, 0d0)
+C-----
+C FLOAT OVERFLOW
+c	CALL ERRSET(72,.TRUE.,.FALSE.,.FALSE.,.FALSE.,15)
+C FLOAT UNDERFLOW
+c	CALL ERRSET(74,.TRUE.,.FALSE.,.FALSE.,.FALSE.,15)
+C EXP TOO SMALL
+c	CALL ERRSET(89,.TRUE.,.FALSE.,.FALSE.,.FALSE.,15)
+C EXP TOO LRGE
+c	CALL ERRSET(88,.TRUE.,.FALSE.,.FALSE.,.FALSE.,15)
+C-----
+      LERROR=1.E30
+	IFLAG=0
+C     COEFFICIENTS
+C-----
+C ASSIGN EXTERNAL PARAMETERS
+	M=IM
+	NN=INN
+	N=IN
+	NT=INTA
+	EPSI=XEPSI
+	DO 8100 I=1,M
+		V(I)=XV(I)
+		YD(I,1)=XYD(I)
+8100	CONTINUE
+C-----
+      T=1.0
+      A=1.0
+      B=0.5
+      G=2.0
+      ICOUNT=0
+      INDEX=1
+      IQ=3*N
+      DO 140 J=1,N
+140   X(1,J)=1.0
+160   DO 172 J=1,N
+172   XF(J)=X(1,J)
+      CALL FVAL
+      Y(1)=YVAL
+      SOLD=YVAL
+C---- CONSTRUCT SIMPLEX
+      EN=N
+      PN=(SQRT(EN+1.0)-1.0+EN)/(EN*SQRT(2.0))*T
+      QN=(SQRT(EN+1.0)-1.0)/(EN*SQRT(2.0))*T
+      NP1=N+1
+      DO 305 I=2,NP1
+      INDEX=I
+      DO 300 J=1,N
+      EJ=0.0
+      EI=0.0
+      IF(I-1.NE.J) EJ=1.0
+      IF(I-1.EQ.J) EI=1.0
+      X(I,J)=X(1,J)+EI*PN+EJ*QN
+300   XF(J)=X(I,J)
+      CALL FVAL
+305   Y(I)=YVAL
+C---- DETERMINE MAX XH
+310   IH=1
+      DO 350 J=1,NP1
+      IF(Y(IH).GE.Y(J)) GOTO 350
+      IH=J
+350   CONTINUE
+C---- DETERMINE SECOND MAX XS
+      IS=1
+      IF(IH.NE.1) GOTO 470
+      IS=2
+470   CONTINUE
+      DO 420 J=1,NP1
+      IF(J.EQ.IH) GOTO 420
+      IF(Y(IS).GE.Y(J)) GOTO 420
+      IS=J
+420   CONTINUE
+C---- DETERMINE MIN XL
+      IL=1
+      DO 480 J=1,NP1
+      IF(Y(IL).LE.Y(J)) GOTO 480
+      IL=J
+480   CONTINUE
+C---- COMPUTE CENTROID
+      DO 510 J=1,N
+510   SUMC(J)=0.0
+      EN=N
+      DO 570 J=1,N
+      DO 560 I=1,NP1
+      IF(I.EQ.IH) GOTO 560
+      SUMC(J)=SUMC(J)+X(I,J)
+560   CONTINUE
+570   XC(J)=1.0/EN*SUMC(J)
+      DO 573 J=1,N
+573   XF(J)=XC(J)
+      CALL FVAL
+      YBAR=YVAL
+      SUM=0.0
+      DO 577 I=1,NP1
+  577 SUM = SUM + ((Y(I)-YBAR)/YBAR)**2
+      ICOUNT=ICOUNT+1
+      ERROR=SQRT(SUM/EN)
+      IQ=IQ-1
+      IF(IQ.EQ.-1) CALL QADFIT
+      IF(IFLAG.EQ.1) GOTO 1990
+      IF(ERROR.LE.EPSI) GOTO 1990
+      IF(ABS(LERROR-ERROR).LT.EPSI) GO TO 1990
+      LERROR=ERROR
+C---- DO A REFLECTION
+      DO 600 J=1,N
+600   XR(J)=(1.0+A)*XC(J)-A*X(IH,J)
+      DO 610 J=1,N
+610   XF(J)=XR(J)
+      INDEX=N+2
+      CALL FVAL
+      YXR=YVAL
+      IF(YXR.GE.Y(IL)) GOTO 750
+C---- DO A EXPANSION
+      DO 660 J=1,N
+660   XE(J)=G*XR(J)+(1.0-G)*XC(J)
+      DO 680 J=1,N
+680   XF(J)=XE(J)
+      INDEX=N+3
+      CALL FVAL
+      YXE=YVAL
+      IF(YXE.GT.Y(IL)) GOTO 760
+      DO 730 J=1,N
+730   X(IH,J)=XE(J)
+      Y(IH)=YXE
+      NP3=N+3
+      DO 735 K=1,M
+735   F(IH,K)=F(NP3,K)
+      GOTO 310
+750   IF(YXR.GT.Y(IS)) GOTO 800
+760   DO 780 J=1,N
+780   X(IH,J)=XR(J)
+      Y(IH)=YXR
+      NP2=N+2
+      DO 785 K=1,M
+785   F(IH,K)=F(NP2,K)
+      GOTO 310
+800   IF(YXR.GT.Y(IH)) GOTO 830
+      DO 820 J=1,N
+820   X(IH,J)=XR(J)
+C---- DO A CONTRACTION
+830   DO 840 J=1,N
+840   XCO(J)=B*X(IH,J)+(1.0-B)*XC(J)
+      DO 860 J=1,N
+860   XF(J)=XCO(J)
+      INDEX=N+2
+      CALL FVAL
+      YXCO=YVAL
+      IF(YXCO.GT.Y(IH)) GOTO 930
+      DO 910 J=1,N
+910   X(IH,J)=XCO(J)
+      Y(IH)=YXCO
+      NP2=N+2
+      DO 915 K=1,M
+915   F(IH,K)=F(NP2,K)
+      GOTO 310
+930   DO 960 I=1,NP1
+      INDEX=I
+      DO 955 J=1,N
+950   X(I,J)=0.5*(X(I,J)+X(IL,J))
+955   XF(J)=X(I,J)
+      CALL FVAL
+960   Y(I)=YVAL
+C---- HAS A MIN BEEN REACHED?
+      GOTO 310
+1990  DO 1594 J=1,N
+	PARS(J)=X(IL,J)
+1594  XF(J)=X(IL,J)
+      CALL FVAL
+	DO 1595 I=1,NT
+1595	BPARS(I)=BB(I)
+      CALL INDEXD
+	RETURN
+      END
+C---------------------------------------------------------------------
+      SUBROUTINE MATIN
+C---- DETERMINE INVERSE OF MATRIX
+      COMMON /NLC/ EPSI,IFLAG,IL,IQ,INDEX,F(15,120),M,N
+      COMMON /NLC/ SOLD,Y(20),YVAL,XF(11),X(11,11),V(120),YD(120,10)
+      COMMON /NLC/ GG(11,11),GINV(11,11),EM(15,120),BB(5),NT
+      COMMON /NLC/ GA(120,5),NN
+      DIMENSION E(15,120),EN(20),T(20),Z(11,11),YY(20)
+      EQUIVALENCE (EM(1,1),E(1,1))
+      DO 20 I=1,N
+      DO 20 J=1,N
+      IF(I.EQ.J) GOTO 10
+      Z(I,J)=0.0
+      GOTO 20
+10    Z(I,J)=1.0
+20    CONTINUE
+      DO 120 J0=1,N
+      I0=J0
+      DO 30 I=1,N
+30    YY(I)=GG(I,J0)
+      DO 40 I=1,N
+      EN(I) = 0.
+      T(I)=0.0
+      DO 40 J=1,N
+40    T(I)=T(I)+Z(I,J)*YY(J)
+      IF(T(J0).EQ.0.) GO TO 65
+      DO 60 J=1,N
+      IF(J.EQ.J0) GOTO 50
+      EN(J)=-T(J)/T(J0)
+      GOTO 60
+50    EN(J)=1./T(J0)
+60    CONTINUE
+   65 DO 80 I = 1,N
+      DO 80 J=1,N
+      IF (I.EQ.J) GOTO 70
+      E(I,J)=0.0
+      GOTO 80
+70    E(I,J)=1.0
+80    CONTINUE
+      DO 90 J=1,N
+90    E(J,J0)=EN(J)
+      DO 100 K=1,N
+      DO 100 I=1,N
+      GINV(K,I)=0.0
+      DO 100 J=1,N
+100   GINV(K,I)=GINV(K,I)+E(K,J)*Z(J,I)
+      DO 110 J=1,N
+      DO 110 I=1,N
+110   Z(I,J)=GINV(I,J)
+120   CONTINUE
+      RETURN
+      END
+C-------------------------------------------------------------------------
+      SUBROUTINE QADFIT
+      COMMON /NLC/ EPSI,IFLAG,IL,IQ,INDEX,F(15,120),M,N
+      COMMON /NLC/ SOLD,Y(20),YVAL,XF(11),X(11,11),V(120),YD(120,10)
+      COMMON /NLC/ GG(11,11),GINV(11,11),EM(15,120),BB(5),NT
+      COMMON /NLC/ GA(120,5),NN
+      DIMENSION A(11,11),DELX(20),E(20),F0(20)
+      NP1=N+1
+C---- QUADRATIC COEFFICIENTS
+      II=0
+      DO 30 K=1,M
+      II=0
+      DO 30 I=1,NP1
+      IF(I.EQ.IL) GOTO 30
+      II=II+1
+      EM(II,K)=F(I,K)-F(IL,K)
+30    CONTINUE
+      DO 50 I=1,N
+      F0(I)=0.0
+      DO 50 K=1,M
+50    F0(I)=F0(I)-F(IL,K)*EM(I,K)
+C---- ELEMENTS OF THE MATRIX GAMMA,G
+      DO 70 I=1,N
+      DO 70 J=1,N
+      GG(I,J)=0.0
+      DO 70 K=1,M
+70    GG(I,J)=GG(I,J)+EM(I,K)*EM(J,K)
+      CALL MATIN
+      DO 80 I=1,N
+      E(I)=0.0
+      DO 80 J=1,N
+80    E(I)=E(I)+GINV(I,J)*F0(J)
+C---- DEFINE THE SCALING MATRIX A
+      II=0
+      DO 101 I=1,NP1
+      IF(I.EQ.IL) GOTO 101
+      II=II+1
+      DO 100 J=1,N
+      A(II,J)=X(I,J)-X(IL,J)
+100   CONTINUE
+101   CONTINUE
+C---- DETERMINE DEL X
+      DO 110 I=1,N
+      DELX(I)=0.0
+      DO 110 J=1,N
+110   DELX(I)=DELX(I)+A(J,I)*E(J)
+      DO 120 J=1,N
+120   XF(J)=X(IL,J)+DELX(J)
+      INDEX=N+2
+      CALL FVAL
+      IF(Y(IL).LT.YVAL) GOTO 140
+      TEMP=ABS(1-SOLD/YVAL)
+      IF(TEMP.EQ.1) GOTO 150
+      IF(TEMP.LE.EPSI) GOTO 150
+      SOLD=YVAL
+      DO 130 J=1,N
+130   X(IL,J)=XF(J)
+      NP2=N+2
+      DO 135 K=1,M
+135   F(IL,K)=F(NP2,K)
+      IFLAG=2
+      IQ=(3*N)/2
+      GOTO 160
+140   IFLAG=2
+      IQ=3*N
+      GOTO 160
+150   IFLAG=1
+      DO 155 J=1,N
+155   X(IL,J)=XF(J)
+      Y(IL)=YVAL
+160   RETURN
+      END
+C----------------------------------------------------------------------
+      SUBROUTINE  INDEXD
+      COMMON /NLC/ EPSI,IFLAG,IL,IQ,INDEX,F(15,120),M,N
+      COMMON /NLC/ SOLD,Y(20),YVAL,XF(11),X(11,11),V(120),YD(120,10)
+      COMMON /NLC/ GG(11,11),GINV(11,11),EM(15,120),BB(5),NT
+      COMMON /NLC/ GA(120,5),NN
+	COMMON /NLCOUT/ FF(120),PARS(10),BPARS(10)
+      SUM=0.0
+      DO 200 I=1,M
+200   SUM=SUM+V(I)
+      XM=M
+      YBAR=SUM/XM
+      SST=0.0
+      DO 240 I=1,M
+240   SST=SST+(V(I)-YBAR)**2
+      SSR=0.0
+      DO 280 I=1,M
+      FF(I)=0.0
+      DO 260  J=1,NT
+260   FF(I)=BB(J)*GA(I,J)+FF(I)
+280   SSR=SSR+(FF(I)-V(I))**2
+      XINDX=1-SSR/SST
+      SIGMAR=SQRT(SSR/XM)
+      DO 300 I=1,M
+      DIFF=FF(I)-V(I)
+      IF(V(I).EQ.0.) GO TO 295
+      DIFF = DIFF*100./V(I)
+      GO TO 300
+295   DIFF=0.
+300   CONTINUE
+C
+C----  WRITE(1) (FF(I),I=1,M)
+C----  WRITE(1) (V(I),I=1,M)
+      RETURN
+      END
+C---------------------------------------------------------------------------
+      SUBROUTINE FVAL
+      COMMON /NLC/ EPSI,IFLAG,IL,IQ,INDEX,F(15,120),M,N
+      COMMON /NLC/ SOLD,Y(20),YVAL,XF(11),X(11,11),V(120),YD(120,10)
+      COMMON /NLC/ GG(11,11),GINV(11,11),EM(15,120),BB(5),NT
+      COMMON /NLC/ GA(120,5),NN
+      DIMENSION GTGA(11,11),GT(5,120),GGG(5,120),B(5)
+      DIMENSION G(120,5),A(11),TR(5),XP(11)
+      EQUIVALENCE (GG(1,1),GTGA(1,1)),(BB(1),B(1)),(XF(1),A(1)),
+     *(G(1,1),GA(1,1))
+      DO 200 I=1,M
+      DO 100 J=1,NN
+100   XP(J)=YD(I,J)
+C
+C---- LOCATION OF TRANSFORMS
+      CALL TRANS(TR,A,XP)
+C
+      DO 110 J=1,NT
+110   GA(I,J)=TR(J)
+200   CONTINUE
+      DO 230 J=1,NT
+      DO 230 I=1,M
+230   GT(J,I)=GA(I,J)
+      DO 280 K=1,NT
+      DO 280 I=1,NT
+      GTGA(K,I)=0.0
+      DO 280 J=1,M
+280   GTGA(K,I)=GTGA(K,I)+GT(K,J)*GA(J,I)
+      HOLD=N
+      N=NT
+      CALL MATIN
+      N=HOLD
+      DO 350 K=1,NT
+      DO 350 I=1,M
+      GGG(K,I)=0.0
+      DO 350 J=1,NT
+350   GGG(K,I)=GGG(K,I)+GINV(K,J)*GT(J,I)
+      DO 400 K=1,NT
+      B(K)=0.0
+      DO 400 J=1,M
+400   B(K)=B(K)+GGG(K,J)*V(J)
+      YVAL=0.0
+      DO 460 I=1,M
+      FF=0.0
+      DO 240 J=1,NT
+240   FF=B(J)*GA(I,J)+FF
+      F(INDEX,I)=V(I)-FF
+460   YVAL=(V(I)-FF)**2+YVAL
+      RETURN
+      END