From 40e5a5811c6ffce9b0974e93cdd927cbcf60c157 Mon Sep 17 00:00:00 2001 From: Joe Hunkeler Date: Tue, 11 Aug 2015 16:51:37 -0400 Subject: Repatch (from linux) of OSX IRAF --- noao/onedspec/fortran/nlcfit.f | 400 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 400 insertions(+) create mode 100644 noao/onedspec/fortran/nlcfit.f (limited to 'noao/onedspec/fortran/nlcfit.f') 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 -- cgit