C @(#)necripp.for	17.1.1.1 (ESO-DMD) 01/25/02 17:51:31
C===========================================================================
C Copyright (C) 1995 European Southern Observatory (ESO)
C
C This program is free software; you can redistribute it and/or 
C modify it under the terms of the GNU General Public License as 
C published by the Free Software Foundation; either version 2 of 
C the License, or (at your option) any later version.
C
C This program is distributed in the hope that it will be useful,
C but WITHOUT ANY WARRANTY; without even the implied warranty of
C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
C GNU General Public License for more details.
C
C You should have received a copy of the GNU General Public 
C License along with this program; if not, write to the Free 
C Software Foundation, Inc., 675 Massachusetts Ave, Cambridge, 
C MA 02139, USA.
C
C Correspondence concerning ESO-MIDAS should be addressed as follows:
C	Internet e-mail: midas@eso.org
C	Postal address: European Southern Observatory
C			Data Management Division 
C			Karl-Schwarzschild-Strasse 2
C			D 85748 Garching bei Muenchen 
C			GERMANY
C===========================================================================
C
      PROGRAM ECHRIP
C
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C 
C.COPYRIGHT: Copyright (c) 1987 European Southern Observatory,
C                                         all rights reserved
C
C.VERSION: 1.0  ESO-FORTRAN Conversion, AA  19:01 - 6 DEC 1987
C
C.LANGUAGE: F77+ESOext
C
C.AUTHOR: D.PONZ
C
C.IDENTIFICATION
C
C  ECHRIPP  version 1.5 840820
C           version 1.6 901009  P.Ballester   Option to freeze K,Alpha
C
C.KEYWORDS
C
C  ECHELLE, CASPEC, BLAZE FUNCTION
C
C.PURPOSE
C
C  COMPUTE THE ECHELLE CONSTANTS ALPHA AND K FOR EACH ECHELLE ORDER
C  (only if the flag is not set to FREEZE)
C
C.ALGORITHM
C
C  IN DEVELOPMENT ...
C
C.INPUT/OUTPUT
C
C   RIPPL/ECHELLE INPUT OUTPUT  LAMBDA_MIN,DELTA_LAMBDA,K0,A0 
C                   -- in overlapping region
C   RIPPL/ECHELLE INPUT OUTPUT  K0,A0   -- fit SINC SQ.
C
C.VERSION
C 
C 010706	last modif 
C 
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C 
      IMPLICIT NONE
C 
      INTEGER        MADRID
      INTEGER        I,ISTAT
      INTEGER        METHOD,NAXISA,OUTIMA,KUN,KNUL,INIMA
      INTEGER        NPIXA(2),NPTOT(100),NORDER(100)
C 
      INTEGER*8      IPNTRA,IPNTRB
C 
      DOUBLE PRECISION  WSTR(100),STEPA(2),STARTA(2)
      DOUBLE PRECISION  WDELTA,WINIT,K0,ALPHA0
C 
      REAL              CUT(4),XMIN,XMAX,PARS(4)
C 
      CHARACTER*8       TABLE
      CHARACTER*60      INFRA,OUTFRE
      CHARACTER         IDENTA*72,CUNITA*64,FREEZE*12
C
      INCLUDE 'MID_INCLUDE:ST_DEF.INC'
      COMMON /VMR/MADRID(1)
      INCLUDE 'MID_INCLUDE:ST_DAT.INC'
C
C ... INITIALIZE SYSTEM
C
      CALL STSPRO('ECHRIPP')
      CALL STKRDC('P1',1,1,60,I,INFRA,KUN,KNUL,ISTAT)
      CALL STKRDC('P2',1,1,60,I,OUTFRE,KUN,KNUL,ISTAT)
      CALL CLNFRA(INFRA,INFRA,0)
      CALL CLNFRA(OUTFRE,OUTFRE,0)
      CALL STKRDR('INPUTR',1,4,I,PARS,KUN,KNUL,ISTAT)
      TABLE = 'riptab'
      CALL STKRDC ('P5',1,1,12,I,FREEZE,KUN,KNUL,ISTAT)
      IF (PARS(3).LT.1.0) THEN
          K0     = PARS(1)
          ALPHA0 = PARS(2)
          METHOD = 0
          CALL STTPUT(' Ripple correction. Method: SINC',ISTAT)          
      ELSE
          WINIT  = PARS(1)
          WDELTA = PARS(2)
          K0     = PARS(3)
          ALPHA0 = PARS(4)
          METHOD = 1
          CALL STTPUT(' Ripple correction. Method: OVER',ISTAT)          
      END IF
C
C ... MAP INPUT/OUTPUT FRAMES
C
      CALL STIGET(INFRA,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,
     .            3,NAXISA,NPIXA,STARTA,STEPA,IDENTA,CUNITA,
     .            IPNTRA,INIMA,ISTAT)
      CALL STDRDR(INIMA,'LHCUTS',1,4,I,CUT,KUN,KNUL,ISTAT)
      CALL STDRDD(INIMA,'WSTART',1,NPIXA(2),I,WSTR,
     .            KUN,KNUL,ISTAT)
      CALL STDRDI(INIMA,'NPTOT',1,NPIXA(2),I,NPTOT,
     .            KUN,KNUL,ISTAT)
      CALL STDRDI(INIMA,'NORDER',1,NPIXA(2),I,NORDER,
     .            KUN,KNUL,ISTAT)
      CALL STIPUT(OUTFRE,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,
     .            NAXISA,NPIXA,STARTA,STEPA,IDENTA,
     +            CUNITA,IPNTRB,OUTIMA,ISTAT)
      IF (METHOD.EQ.0) THEN
          CALL RIPPL0(NPIXA(1),NPIXA(2),MADRID(IPNTRA),MADRID(IPNTRB),
     +                 WSTR,STEPA(1),NPTOT,NORDER,K0,ALPHA0,XMIN,XMAX,
     +                 TABLE,FREEZE)

      ELSE
          CALL RIPPL1(NPIXA(1),NPIXA(2),MADRID(IPNTRA),MADRID(IPNTRB),
     +                 WSTR,STEPA(1),NPTOT,NORDER,WINIT,WDELTA,K0,
     +                 ALPHA0,XMIN,XMAX,TABLE,FREEZE)
          XMIN   = CUT(3)
          XMIN   = CUT(4)
      END IF
C
C ... WRITE PROCESS DESCRIPTORS
C
      CUT(1) = XMIN
      CUT(2) = XMAX
      CUT(3) = XMIN
      CUT(4) = XMAX
      CALL DSCUPT(INIMA,OUTIMA,' ',ISTAT)
      CALL STDWRR(OUTIMA,'LHCUTS',CUT,1,4,KUN,ISTAT)
      CALL STSEPI
      END

      SUBROUTINE RIPPL1(NX,NY,X,Y,WS,WSTEP,NP,NO,W1,DW,K0,A0,XMIN,XMAX,
     +                   TABLE,FZ)
C
C CORRECT FOR BLAZE EFFECT USING THE OVERLAPPING REGION
C
C INPUT ARGUMENTS
C NX MAXIMUM NUMBER OF SAMPLES PER ORDER
C NY NUMBER OF ORDERS
C X  INPUT IMAGE
C WS STARTING WAVELENGTH PER EACH ORDER
C WSTEP WAVELENGTH STEP
C NP NUMBER OF USEFUL SAMPLES PER ORDER
C NO SPECTRAL ORDER NUMBER
C W1 STARTING OFFSET PER ORDER TO DEFINE THE OVERLAPPING REGION
C DW INTERVAL OF THE OVERLAPPING REGION
C K0 INITIAL GUESS FOR THE CONSTANT K
C A0 INITIAL GUESS FOR ALPHA
C FZ FLAG TO FREEZE THE VALUES OF K and ALPHA
C
C OUTPUT ARGUMENTS
C Y    CORRECTED IMAGE
C XMIN MINIMUN AND
C XMAX MAXIMUM VALUES
C
      IMPLICIT NONE
C
      INTEGER          NX, NY
      INTEGER          NP(NY),NO(NY),IERR(2),ICOL(7), ISEQ
      INTEGER          I, TID, ISTAT, IORDER, I1, ILR, I2
      INTEGER          NN, NK
      REAL             X(NX,NY), Y(NX,NY)
      REAL             VAL(7), XMIN, XMAX
      DOUBLE PRECISION WS(NY), WSTEP, K0, A0
      DOUBLE PRECISION W1, DW, K(2), KM, A(2), AM
      CHARACTER*16     LABEL(7),UNIT(7),FORM
      CHARACTER*12     FZ
      CHARACTER*80     LINE
      CHARACTER*(*)    TABLE
C
      INCLUDE 'MID_INCLUDE:ST_DEF.INC'
      INCLUDE 'MID_INCLUDE:ST_DAT.INC'
C
      DATA FORM/'G12.4'/
      DATA LABEL(1)/'ORDER'/,UNIT(1)/'UNITLESS'/
      DATA LABEL(2)/'K'/,UNIT(2)/'UNITLESS'/
      DATA LABEL(3)/'ALPHA'/,UNIT(3)/'UNITLESS'/
      DATA LABEL(4)/'KFIT'/,UNIT(4)/'UNITLESS'/
      DATA LABEL(5)/'AFIT'/,UNIT(5)/'UNITLESS'/
      DATA LABEL(6)/'YMAX'/,UNIT(6)/'FLUX'/
      DATA LABEL(7)/'STATUS'/,UNIT(7)/'UNITLESS'/

C
      ISEQ   = 0
      XMIN   = 1.E30
      XMAX   = -XMIN
C
C ... INITIALIZE TEMPORARY OUTPUT
C
      CALL TBTINI(TABLE,F_TRANS,F_O_MODE,8,NY,TID,ISTAT)
      DO 10 I = 1,6
          CALL TBCINI(TID,D_R4_FORMAT,1,FORM,UNIT(I),
     .                LABEL(I),ICOL(I),ISTAT)
   10 CONTINUE
C
C ... ITERATION ON ALL CONSECUTIVE PAIRS OF ORDERS
C
      CALL STTPUT(' ORDER NUMER  FITTED K ',ISTAT)
      CALL STTPUT(' ----------- ----------',ISTAT)
      DO 30 IORDER = 2,NY - 1

      IF (FZ(1:1).EQ.'F'.OR.FZ(1:1).EQ.'f') THEN

          KM = K0
          AM = A0

      ELSE

          I1     = IORDER - 2
C
C ... CONSIDER LEFT AND RIGHT OVERLAPPING REGIONS
C
          DO 20 ILR = 0,1
              I1     = I1 + 1
              I2     = I1 + 1
              CALL DELTAK(WS(I1),WSTEP,NP(I1),X(1,I1),NO(I1),WS(I2),
     +                    NP(I2),X(1,I2),NO(I2),W1,DW,K0,A0,K(ILR+1),
     +                    A(ILR+1),IERR(ILR+1))
   20     CONTINUE
C
C ... COMPUTE MEAN VALUE
C
          KM     = 0.5* (K(1)+K(2))
          AM     = 0.5* (A(1)+A(2))

      ENDIF   ! End of the test on FZ

C   WRITE(*,*) IORDER,KM,AM
C
C ... CORRECT FOR RIPPL
C
          CALL RIPCOR(WS(IORDER),WSTEP,NP(IORDER),X(1,IORDER),
     +                Y(1,IORDER),NO(IORDER),KM,AM)
          NN = NP(IORDER)
          DO 35 NK = 1, NN
             IF (XMAX.LT.Y(NK,IORDER)) XMAX = Y(NK,IORDER)
             IF (XMIN.GT.Y(NK,IORDER)) XMIN = Y(NK,IORDER)
35        CONTINUE
          IF (IORDER.EQ.2) THEN
              CALL RIPCOR(WS(1),WSTEP,NP(1),X(1,1),Y(1,1),NO(1),KM,AM)
              VAL(1) = NO(1)
              VAL(2) = KM
              VAL(3) = AM
              VAL(4) = KM
              VAL(5) = AM
              ISEQ   = ISEQ + 1
              CALL TBRWRR(TID,ISEQ,6,ICOL,VAL,ISTAT)
          END IF
          VAL(1) = NO(IORDER)
          VAL(2) = KM
          VAL(3) = AM
          VAL(4) = KM
          VAL(5) = AM
          ISEQ   = ISEQ + 1
          CALL TBRWRR(TID,ISEQ,6,ICOL,VAL,ISTAT)
          IF (IORDER.EQ.NY-1) THEN
              CALL RIPCOR(WS(NY),WSTEP,NP(NY),X(1,NY),Y(1,NY),NO(NY),KM,
     +                    AM)
              VAL(1) = NO(NY)
              VAL(2) = KM
              VAL(3) = AM
              VAL(4) = KM
              VAL(5) = AM
              ISEQ   = ISEQ + 1
              CALL TBRWRR(TID,ISEQ,6,ICOL,VAL,ISTAT)
          END IF
C
C ... WRITE OUT RESULTS
C
          WRITE (LINE,9000) NO(IORDER),KM
          CALL STTPUT(LINE,ISTAT)
   30 CONTINUE
      CALL TBSINI(TID,ISTAT)
      CALL TBTCLO(TID,ISTAT)
      RETURN
 9000 FORMAT (4X,I8,F12.2)
      END

      SUBROUTINE RIPPL0(NX,NY,X,Y,WS,WSTEP,NP,NO,K0,A0,XMIN,XMAX,
     +                                                  TABLE,FZ)
C
C CORRECT FOR BLAZE EFFECT BY FITTING SINC SQ. FUNCTION
C
C INPUT ARGUMENTS
C NX MAXIMUM NUMBER OF SAMPLES PER ORDER
C NY NUMBER OF ORDERS
C X  INPUT IMAGE
C WS STARTING WAVELENGTH PER EACH ORDER
C WSTEP WAVELENGTH STEP
C NP NUMBER OF USEFUL SAMPLES PER ORDER
C NO SPECTRAL ORDER NUMBER
C K0 INITIAL GUESS FOR THE CONSTANT K
C A0 INITIAL GUESS FOR ALPHA
C FZ FLAG TO FREEZE THE VALUES OF K AND ALPHA
C
C OUTPUT ARGUMENTS
C Y    CORRECTED IMAGE
C XMIN MINIMUN AND
C XMAX MAXIMUM VALUES
C
      IMPLICIT NONE
C
      INTEGER          NX, NY
      INTEGER          NP(NY),NO(NY),IERR(100),ICOL(7)
      INTEGER          ISEQ, TID, ISTAT, I, IORDER, NNN, I1, N
      INTEGER          NN, NK
      REAL             X(NX,NY), Y(NX,NY), XMIN, XMAX
      REAL             VAL(7), YMAX(100)
      DOUBLE PRECISION WS(NY), WSTEP, K0, A0, W1
      DOUBLE PRECISION A(100), AFIT(100), K(100), KFIT(100)
      CHARACTER*12     FZ
      CHARACTER*16     LABEL(7),UNIT(7),FORM
      CHARACTER*(*)    TABLE
C
C
      INCLUDE 'MID_INCLUDE:ST_DEF.INC'
      INCLUDE 'MID_INCLUDE:ST_DAT.INC'
C
      DATA FORM/'G12.4'/
      DATA LABEL(1)/'ORDER'/,UNIT(1)/'UNITLESS'/
      DATA LABEL(2)/'K'/,UNIT(2)/'UNITLESS'/
      DATA LABEL(3)/'ALPHA'/,UNIT(3)/'UNITLESS'/
      DATA LABEL(4)/'KFIT'/,UNIT(4)/'UNITLESS'/
      DATA LABEL(5)/'AFIT'/,UNIT(5)/'UNITLESS'/
      DATA LABEL(6)/'YMAX'/,UNIT(6)/'FLUX'/
      DATA LABEL(7)/'STATUS'/,UNIT(7)/'UNITLESS'/
C
      ISEQ   = 0
      XMIN   = 1.E30
      XMAX   = -XMIN
C
C ... INITIALIZE TEMPORARY OUTPUT
C
      CALL TBTINI(TABLE,F_TRANS,F_O_MODE,8,NY,TID,ISTAT)
      DO 10 I = 1,6
          CALL TBCINI(TID,D_R4_FORMAT,1,FORM,UNIT(I),
     .                LABEL(I),ICOL(I),ISTAT)
   10 CONTINUE
C
C ... ITERATION ON ORDERS
C
      DO 20 IORDER = 1,NY
          NNN    = 10./WSTEP
          I1     = 1 + NNN/2
          W1     = WS(IORDER) + (I1-1)*WSTEP
          N      = NP(IORDER) - NNN

          CALL FITPAR(W1,WSTEP,N,X(I1,IORDER),NO(IORDER),K0,A0,
     +                K(IORDER),A(IORDER),YMAX(IORDER),IERR(IORDER))

          IF (FZ(1:1).EQ.'F'.OR.FZ(1:1).EQ.'f') THEN
             K(IORDER) = K0
             A(IORDER) = A0
          ENDIF

   20 CONTINUE
C
C ... COMPUTE FITTED RESULTS
C
      CALL COMFIT(NY,NO,K,IERR,KFIT,1)
      CALL COMFIT(NY,NO,A,IERR,AFIT,0)
C
C ... CORRECT FOR BLAZE FUNCTION
C
      DO 30 IORDER = 1,NY
          CALL RIPCOR(WS(IORDER),WSTEP,NP(IORDER),X(1,IORDER),
     +                Y(1,IORDER),NO(IORDER),KFIT(IORDER),AFIT(IORDER))
          VAL(1) = NO(IORDER)
          VAL(2) = K(IORDER)
          VAL(3) = A(IORDER)
          VAL(4) = KFIT(IORDER)
          VAL(5) = AFIT(IORDER)
          VAL(6) = YMAX(IORDER)
          VAL(7) = IERR(IORDER)
          ISEQ   = ISEQ + 1
          CALL TBRWRR(TID,ISEQ,6,ICOL,VAL,ISTAT)
          NN = NP(IORDER)
          DO 35 NK = 1, NN
             IF (XMAX.LT.Y(NK,IORDER)) XMAX = Y(NK,IORDER)
             IF (XMIN.GT.Y(NK,IORDER)) XMIN = Y(NK,IORDER)
35        CONTINUE
C
C ... WRITE OUT RESULTS
C
   30 CONTINUE
      CALL TBSINI(TID,ISTAT)
      CALL TBTCLO(TID,ISTAT)
      RETURN
      END

      SUBROUTINE RIPCOR(W1,WS,NP,X,Y,M,K,A)
C
C CORRECT FOR RIPPL FUNCTION
C
C INPUT ARGUMENTS
C W1 INITIAL WAVELENGTH
C WS WAVELENGTH STEP
C NP NO. OF POINTS
C X INPUT FLUX
C M ORDER NUMBER
C K CONSTANT
C A ALPHA CONSTANT
C
C OUTPUT ARGUMENTS
C Y CORRECTED FLUX
C
      IMPLICIT NONE
C
      INTEGER          NP, M
      REAL             X(NP),Y(NP)
      DOUBLE PRECISION W1, WS, K, A
      DOUBLE PRECISION PI,DK,DA,DW,DM,R,DX,PA,DC
C
      INTEGER          I
      DATA PI/3.1415926535897932D0/
C
      DK     = K
      DA     = A
      DM     = M
      PA     = PI*A
      DC     = DM/DK
      DO 10 I = 1,NP
          DW     = W1 + (I-1)*WS
          DX     = (PA*DM*DC)* (DW-1.D0/DC)
          IF (DABS(DX).LT.1.E-10) THEN
              Y(I)   = X(I)
          ELSE
              R      = DSIN(DX)/DX
              Y(I)   = X(I)/ (R*R)
          END IF
   10 CONTINUE
      RETURN
      END

      SUBROUTINE DELTAK(WSTR1,WSTEP,N1,X1,NO1,WSTR2,N2,X2,NO2,W1,DW,
     +                  K0,A0,K,A,IERR)
C
C COMPUTE CONTANT K
C
C INPUT ARGUMENTS
C WSTEP REAL*4 WAVELENGTH STEP
C WSTR1 REAL*4 STARTING WAVELENGTH ORDER M+1
C N1 INTG*4 NUMBER OF POINTS IN ORDER M+1
C X1 REAL*4 FLUX IN ORDER M+1
C NO1 INTG*4 ORDER M+1
C WSTR2 REAL*4 STARTING WAVELENGTH ORDER M
C N2 INTG*4 NUMBER OF POINTS IN ORDER M
C X2 REAL*4 FLUX IN ORDER M
C NO2 INTG*4 ORDER M
C W1 OVERLAPPING OFFSET
C DW OVERLAPPING INTERVAL
C K0 INITIAL K VALUE
C A0 ALPHA
C
C OUTPUT ARGUMENTS
C K FINAL K VALUE
C A FINAL A VALUE
C IERR CONVERGENCE FLAG 0 - GOOD VALUE
C      OTHERWISE ITT DOES NOT CONVERGE
C
C
      IMPLICIT NONE
C
      INTEGER          N1, N2, NO1, NO2, NPIX, I1, I2, I, IERR
      INTEGER          NCOUNT
      REAL             X1(N1), X2(N2)
      DOUBLE PRECISION K0, K, A, A0
      DOUBLE PRECISION WSTR1, WSTEP, WSTR2, W1, DW
      DOUBLE PRECISION KVAL,AVERP,AVERM,WL,XP,RP,XM,RM
      DOUBLE PRECISION FRAT,R1,R2,PA,PI
      DOUBLE PRECISION XX1,XX2,S1,S2,C1,C2,DK
      DOUBLE PRECISION DM1,DM2,C11,C12,C21,C22
      DOUBLE PRECISION AVAL
      DATA PI/3.1415926535897932D0/
C
      NPIX   = DW/WSTEP
      I1     = (WSTR2+W1-WSTR1)/WSTEP
      I2     = W1/WSTEP
      AVAL   = A0
      PA     = AVAL*PI
      KVAL   = K0
      DM1    = NO1
      DM2    = NO2
      NCOUNT = 0
   10 CONTINUE
      AVERP  = 0.D0
      AVERM  = 0.D0
      C11    = PA*DM1* (DM1/KVAL)
      C12    = PA*DM2* (DM2/KVAL)
      C21    = KVAL/DM1
      C22    = KVAL/DM2
C
C ... COMPUTE RATIO OF THE AVERAGED CORRECTED FLUXES
C
      DO 20 I = 1,NPIX
          WL     = WSTR2 + (I2+I)*WSTEP
          XP     = C12* (WL-C22)
          RP     = DSIN(XP)/XP
          AVERP  = AVERP + X2(I2+I)/ (RP*RP)/NPIX
          XM     = C11* (WL-C21)
          RM     = DSIN(XM)/XM
          AVERM  = AVERM + X1(I1+I)/ (RM*RM)/NPIX
   20 CONTINUE
      FRAT   = AVERP/AVERM
C
C ... COMPUTE INCREMENT IN K
C
      WL     = WSTR2 + W1 + DW*0.5
      XX1    = C11* (WL-C21)
      S1     = DSIN(XX1)
      C1     = DCOS(XX1)
      R1     = 2.D0* (1.D0-XX1*C1/S1)*WL/ (KVAL* (WL-KVAL/NO1))
C T1  = 2.D0*(XX1*C1/S1-1.D0)/AVAL
      XX2    = C12* (WL-C22)
      S2     = DSIN(XX2)
      C2     = DCOS(XX2)
      R2     = 2.D0* (1.D0-XX2*C2/S2)*WL/ (KVAL* (WL-KVAL/NO2))
C T2  = 2.D0*(XX2*C2/S2-1.D0)/AVAL
      DK     = (FRAT-1.0)/ (R2-FRAT*R1)
C DA  = (FRAT-1.0)/(T2-FRAT*T1)
      KVAL   = KVAL + DK
C AVAL   = AVAL + DA
      NCOUNT = NCOUNT + 1
C WRITE(*,*) NCOUNT,KVAL,DK
      IF (DABS(DK).GT.1.D0 .AND. NCOUNT.LT.10) GO TO 10
      IF (NCOUNT.GE.10) THEN
          IERR   = 1
          K      = K0
          A      = A0
      ELSE
          IERR   = 0
          K      = KVAL
          A      = AVAL
      END IF
      RETURN
      END

      SUBROUTINE FITPAR(FSTART,FSTEP,FNPIX,YY,NORDER,K0,A0,K,
     +                  A1,FMAX,IERR)
C
C ITERATE ON THE FITTING PROCESS FOT IMAGE DATA
C THE ITERATION IS FINISHED IF
C - NITER IS REACHED OR
C - CHISQ IS REACHED
C - CONVERGENCE IS ACHIEVED
C
C INPUT PARAMETERS
C YY VALUES OF THE IMAGE
C K0 GUESS FOR K
C A0 INPUT GUESS FOR ALPHA
C
C OUTPUT PARAMETERS
C K  ACTUAL K
C A1 ACTUAL ALPHA
C IERR STATUS RETURN
C
      IMPLICIT NONE
C
      INTEGER          FNPTOT,FNITER,FNPIX
      INTEGER          NITER, NITER1, I, I2, I1, ISTAT, ICOUNT
      INTEGER          J, NACT, IX1, KK, NN, NP, NNP, IERR, NORDER
      REAL             X,YY(FNPIX), FMAX, ACTCH
      REAL             CHISQ, RELAX, XM, CHIPER, CHIOLD
      DOUBLE PRECISION CHI,PP,Q,Y,Y1,Y2, FSTART, FSTEP
      DOUBLE PRECISION AL,B(10),A(10,10), W, K, K0, A1, A0
      DOUBLE PRECISION FDERIV(3),PI
      DOUBLE PRECISION FVALUE(3),FCHISQ,FCCHIS,FRELAX
      CHARACTER*80     LINE
      DATA PI/3.1415926535897932D0/
C
      NITER  = 20
      CHISQ  = 0.
      RELAX  = 0.8
      NITER1 = NITER
      FNPTOT = 3
C
C      WRITE (*,*) ' *** ORDER ',NORDER,' ***'
      XM     = -9999.
      DO 10 I = 1,FNPIX
          XM     = MAX(XM,YY(I))
   10 CONTINUE  
      FVALUE(1) = XM  
      FVALUE(2) = -PI*A0*NORDER  
      FVALUE(3) = -(FVALUE(2)*NORDER)/K0
      I1     = 1
      I2     = FNPIX
      FCCHIS = 9999.
      FCHISQ = CHISQ
      FRELAX = RELAX
      CHI    = 0.D0
      ICOUNT = 0
      CALL STTPUT(' ITERATION  CHI SQ. PERC.CHANGE',ISTAT)
      CALL STTPUT(' ---------  ------- -----------',ISTAT)
      CHIPER = 100.
   15 CONTINUE
          IF ( .NOT. (ICOUNT.LT.NITER1.AND.CHIPER.GT.0.1)) GO TO 100
          DO 30 I = 1,10
              DO 20 J = 1,10
                  A(J,I) = 0.D0
   20         CONTINUE
              B(I)   = 0.D0
   30     CONTINUE
          CHIOLD = CHI
          CHI    = 0.D0
C
C ITERATION ON IMAGE DATA
C
          W      = 1.
          DO 60 IX1 = I1,I2
              X      = FSTART + (IX1-1)*FSTEP
              Y      = YY(IX1)
C
C ...       COMPUTE FUNCTION VALUES AND DERIVATIVES (0 IF FIXED PARAM)
C
              CALL FTFUNC(X,FNPTOT,FVALUE,Y1,FDERIV)
C
C ...       COMPUTE RESIDUALS, WEIGHTED RESIDUALS AND CHI**2
C
              Y2     = Y - Y1
              PP     = DBLE(W)
              Q      = Y2
              CHI    = CHI + Q*Q*PP
C
C ...       BUILD NORMAL EQUATION MATRIX
C
              DO 50 J = 1,FNPTOT
                  AL     = FDERIV(J)
                  DO 40 KK = 1,FNPTOT
                      A(J,KK) = A(J,KK) + AL*PP*FDERIV(KK)
   40             CONTINUE
                  B(J)   = B(J) + PP*Q*AL
   50         CONTINUE
   60     CONTINUE
          NN     = FNPIX
          NP     = FNPTOT
C
C ...     INVERT NORMAL EQUATION MATRIX TO GET THE COVARIANCE MATRIX
C
          NNP    = FNPTOT
          CALL DMATIN(A,FNPTOT,10,ISTAT)
          IF (ISTAT.NE.0) THEN
              CALL STTPUT(' Problems inverting matrix ... ',ISTAT)
              IERR   = 2
              K      = K0
              A1     = A0
              FMAX   = FVALUE(1)
              RETURN
          END IF
C
C ... COMPUTE REDUCED CHI**2 AND CORRECTIONS FOR THE UNKNOWNS
C
          CHI    = DMAX1(CHI/DBLE(NN-NP),0.D0)
          DO 80 I = 1,FNPTOT
              Q      = 0.0D0
              DO 70 J = 1,FNPTOT
                  Q      = Q + A(I,J)*B(J)
   70         CONTINUE
C
C ... COMPUTE ERRORS OF THE CORRECTIONS (I.E. THE UNKNOWNS)
C ... AND APPLY CORRECTIONS (DO THE ADJUSTMENT!)
C
              FVALUE(I) = FVALUE(I) + RELAX*Q
   80     CONTINUE
          ICOUNT = ICOUNT + 1
          FNITER = FNITER + 1
          FCCHIS = SNGL(DSQRT(CHI))
          IF (ICOUNT.GT.1) THEN
              CHIPER = 100.* (CHIOLD-CHI)/CHIOLD
              WRITE (LINE,9000) ICOUNT,FCCHIS,CHIPER
          ELSE
              WRITE (LINE,9000) ICOUNT,FCCHIS
          END IF
          CALL STTPUT(LINE,ISTAT)
      GOTO 15
  100 CONTINUE
      CALL STTPUT(' ',ISTAT)
      NACT   = FNITER
      ACTCH  = FCCHIS
      A1     = -FVALUE(2)/ (PI*NORDER)
      K      = - (FVALUE(2)*NORDER)/FVALUE(3)
      FMAX   = FVALUE(1)
      IERR   = 0
      IF (CHIPER.GT.0.1) IERR   = 1
      RETURN
 9000 FORMAT (I10,1X,G12.4,3X,F7.3)
C 9010 FORMAT (I10,1X,G12.4,1X,G12.4)
      END

      SUBROUTINE DMATIN(A,N,NDIM,ISTAT)
C
C MATRIX INVERSION USING THE GAUSS-JORDAN METHOD
C  MODIFIED VERSION BY F.L. BAUER AND C. REINSCH
C
      IMPLICIT NONE
      INTEGER N, NDIM, ISTAT
      DOUBLE PRECISION A(NDIM,NDIM)
C
      INTEGER K, KK, I, J
      DOUBLE PRECISION P,Q,H(400)

C
      ISTAT  = 0
      DO 40 K = 1,N
          IF (A(1,1).LE.0.0D0) GO TO 70
          KK     = N + 1 - K
          P      = 1.0/A(1,1)
          DO 20 I = 2,N
              Q      = A(I,1)
              H(I)   = Q*P
              IF (I.LE.KK) H(I)   = -H(I)
              DO 10 J = 2,I
                  A(I-1,J-1) = A(I,J) + Q*H(J)
   10         CONTINUE
   20     CONTINUE
          A(N,N) = P
          DO 30 I = 2,N
              A(N,I-1) = H(I)
   30     CONTINUE
   40 CONTINUE
      DO 60 I = 2,N
          DO 50 J = 1,I
              A(J,I) = A(I,J)
   50     CONTINUE
   60 CONTINUE
      RETURN
   70 ISTAT  = 1
      RETURN
      END

      SUBROUTINE FTFUNC(X,NP,PARAM,Y1,DERIV)
C
C FUNCTION DEFINITION
C
C INPUT PARAMS
C IND INTG NO. OF IND. VARS
C X REAL IND.VARS
C NP INTG NO. OF PARAMS
C PARAM DBLE FUNCTION PARAMETERS
C
C OUTPUT PARAMS
C Y1 DBLE COMPUTED FUNCTION VALUE
C DERIV DBLE COMPUTED DERIVATIVE VALUE
C
      IMPLICIT NONE
      INTEGER          NP
      REAL             X
      DOUBLE PRECISION Y1,PARAM(NP),DERIV(NP)
C
C SINC**2
C Y = P(1)*SINC(P(2)+P(3)*X)**2
C
      DOUBLE PRECISION W, W2
C
      W      = PARAM(3)*X + PARAM(2)
      IF (DABS(W).LT.10.E-20) THEN
          Y1     = PARAM(1)
          DERIV(1) = 1.D0
          DERIV(2) = 0.D0
          DERIV(3) = 0.D0
      ELSE
          W2       = DSIN(W)/W
          DERIV(1) = W2*W2
          Y1       = DERIV(1)*PARAM(1)
          DERIV(2) = 2.D0*PARAM(1)*W2* (DCOS(W)-W2)/W
          DERIV(3) = X*DERIV(2)
      END IF
      RETURN
      END

      SUBROUTINE COMFIT(N,IX,Y,IERR,YFIT,IDEG)
C
C COMPUTE LINEAR FIT
C
C INPUT ARGUMENTS
C N INTG*4 NUMBER OF VALUES
C IX INTG*4 ARRAY WITH THE INDEPENDENT VARIABLE
C Y REAL*4 ARRAY WITH THE DEPENDENT VARIABLE
C IERR INTG*4 ERROR FLAG
C IDEG INTG*4 DEGREE 0 OR 1
C
C OUTPUT ARGUMENTS
C YFIT REAL*4 FITTED VALUES
C
      IMPLICIT NONE
      INTEGER          N, I, NN, IDEG
      INTEGER          IX(N),IERR(N)
      DOUBLE PRECISION Y(N),YFIT(N)
      DOUBLE PRECISION SY,SX,SX2,SXY,B0,B1,YF,YRES,YSIG,THRES
C
      SX     = 0.D0
      SY     = 0.D0
      SX2    = 0.D0
      SXY    = 0.D0
      NN     = 0
C
C ... FIRST ITERATION
C
      DO 10 I = 1,N
          IF (IERR(I).EQ.0) THEN
              NN     = NN + 1
              SY     = SY + Y(I)
              IF (IDEG.EQ.1) THEN
                  SXY    = SXY + IX(I)*Y(I)
                  SX     = SX + IX(I)
                  SX2    = SX2 + IX(I)*IX(I)
              END IF
          END IF
   10 CONTINUE
C
C ... COMPUTE FIRST APPROX
C
      IF (IDEG.EQ.1) THEN
          B1     = (NN*SXY-SX*SY)/ (NN*SX2-SX*SX)
      ELSE
          B1     = 0.D0
      END IF
      B0     = SY/NN - B1* (SX/NN)
C
C ... COMPUTE RESIDUALS
C
      SX     = 0.D0
      SX2    = 0.D0
      DO 20 I = 1,N
          YF     = B1*IX(I) + B0
          YFIT(I) = (Y(I)-YF)/YF
          SX     = SX + YFIT(I)
          SX2    = SX2 + YFIT(I)*YFIT(I)
   20 CONTINUE
      YRES   = SX/N
      YSIG   = DSQRT((SX2-N*YRES*YRES)/ (N-1))
C
C ... SECOND ITERATION
C
      SX     = 0.D0
      SY     = 0.D0
      SX2    = 0.D0
      SXY    = 0.D0
      NN     = 0
      THRES  = 2.D0*YSIG
      DO 30 I = 1,N
          IF (DABS(YFIT(I)-YRES).LE.THRES) THEN
              NN     = NN + 1
              SY     = SY + Y(I)
              IF (IDEG.EQ.1) THEN
                  SXY    = SXY + IX(I)*Y(I)
                  SX     = SX + IX(I)
                  SX2    = SX2 + IX(I)*IX(I)
              END IF
          END IF
   30 CONTINUE
C
C ... COMPUTE FIRST APPROX
C
      IF (IDEG.EQ.1) THEN
          B1     = (NN*SXY-SX*SY)/ (NN*SX2-SX*SX)
      ELSE
          B1     = 0.D0
      END IF
      B0     = SY/NN - B1* (SX/NN)
C
C ... COMPUTE FINAL VALUES
C
      DO 40 I = 1,N
          YFIT(I) = B1*IX(I) + B0
   40 CONTINUE
      RETURN
      END