C @(#)necwlcopt.for	19.1 (ES0-DMD) 02/25/03 14:20:26
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 Massachusetss Ave, Cambridge, 
C MA 02139, USA.
C
C Corresponding 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
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++
C.COPYRIGHT   (C) 1992 European Southern Observatory
C.IDENT       wlcoptim.for
C.AUTHOR      Pascal Ballester,  ESO - Garching
C.KEYWORDS    Spectroscopy, Echelle,
C.PURPOSE     Computes optimal parameters for the 
C             wavelength calibration
C.VERSION     1.0    Creation    1992-MAY-26
C-------------------------------------------------------
C
C
C  Program : Wavelength Calibration Least Squares Solutions
C
C  This program has been partly generated using Macsyma to 
C  solve the following least-square problem:
C
C     Xr(i) = cos(a) * X(i) + bin * sin(a) * Y(i)    (1)
C
C  takes into account a binning factor of the ccd and a rotation
C  of the detector.
C
C     F(Xr(i)) = B * Xr(i) + A = (m+p(i))*lambda(i)  (2)
C
C  is the echelle relation applied in the rotated space and approximated
C  by a linear relation. m is the absolute order number of a line
C  and p(i) is relative order number of the other lines.
C
C  Macsyma has been used to find different optima, like:
C    - The absolute order number of the first identification : m
C    - The binning factor    : bin
C    - In method ANGLE, the angle of rotation : a
C    - In method PAIR, the wavelength of the two pairs : lambda(1), lambda(2)
C
C  The different least-square solutions are used to compute the optimal
C  value of a parameter, assuming that all other parameters are correctly
C  provided. This is a cross-check of the entered values. By computing 
C  the new rms after replacement of a parameter by its optimal value,
C  it is possible to provide a diagnosis of what action would benefit
C  the most to solve the problem.
C
C  INPUT:
C   IN_A/C/1/60   Name of line table
C   IN_B/C/1/60   Method (PAIR or ANGLE)
C   INPUTR/R/1/2  Angle of rotation, Binning ratio.
C
C  OUTPUT:
C   OUTPUTI/I/1/1 Order number of the first line
C   OUTPUTR/R/1/2 Optimal rotation angle, Optimal binning factor
C
        IMPLICIT NONE

        INTEGER  MAXDIM
        PARAMETER (MAXDIM=1500)
        DOUBLE PRECISION    IX(MAXDIM), IY(MAXDIM)
        DOUBLE PRECISION    IP(MAXDIM), IL(MAXDIM)
        DOUBLE PRECISION    SX, SY, SX2, SY2, SXY
C        DOUBLE PRECISION    SML2K, SXMLK, SYMLK, SMLK
        DOUBLE PRECISION    SP2L2K, BUFFER
        DOUBLE PRECISION    SXLK, SXPLK, SYLK, SYPLK, SPLK, SLK
        DOUBLE PRECISION    SL2K, SPL2K, PK, LK, XK, YK, RMSLK
        DOUBLE PRECISION    SML2, SP2L2, SXML, SYML, SML
        DOUBLE PRECISION    SXL, SXPL, SYL, SYPL, SPL, SL
        DOUBLE PRECISION    SL2, SPL2, SBIN
        INTEGER             M, N, K, M2OPT, LOOP, ERROR, P
        DOUBLE PRECISION    BIN, PI, MDP, MINRMS(MAXDIM)
        DOUBLE PRECISION    A, B1
        DOUBLE PRECISION    TGA1, TGA2, AOPT
        REAL                OUTANG, RBIN
        DOUBLE PRECISION    L1,L2,X1,X2,XP1,XP2,L1OPT,L2OPT
        DOUBLE PRECISION    SXR, SXR2, SXRML, DET, COEFA, COEFB, RMS
        DOUBLE PRECISION    L1OLD, L2OLD, RMSL1, RMSL2
        CHARACTER           MODE*5, METHOD*10, TEXT*80
        INTEGER             MOLD, ROW
        DOUBLE PRECISION    RMSM, MRMS, ARMS, BRMS, AOLD, BOLD
        DOUBLE PRECISION    LKOPT(MAXDIM), LKO

        CHARACTER           TABLE*60
        INTEGER             IAV, KUN, KNUL, ISTAT, ICOLX, ICOLW
        INTEGER             ICOLO, ICOLN, TID, NCOL1, NROW1, NSC1
        INTEGER             NACOL, NAROW, MADRID
        LOGICAL             SELECT

        INCLUDE 'MID_INCLUDE:ST_DEF.INC'
        COMMON/VMR/MADRID(1)
        INCLUDE 'MID_INCLUDE:ST_DAT.INC'

C
C Initialisations

C
C ... Init input table
C
      CALL STSPRO('WLCLS')

      CALL STKRDC('IN_A',1,1,60,IAV,TABLE,KUN,KNUL,ISTAT)
      CALL TBTOPN(TABLE,F_U_MODE,TID,ISTAT)
C
      CALL TBIGET(TID,NCOL1,NROW1,NSC1,NACOL,NAROW,ISTAT)
      CALL TBCSER(TID,':X '    ,ICOLX,ISTAT)
      CALL TBCSER(TID,':IDENT' ,ICOLW,ISTAT)
      CALL TBCSER(TID,':ORDER' ,ICOLO,ISTAT)
      CALL TBCSER(TID,':YNEW'  ,ICOLN,ISTAT)

      CALL STKRDC('IN_B',1,1,10,IAV,METHOD,KUN,KNUL,ISTAT)
      CALL STKRDR('INPUTR',1,1,IAV,OUTANG,KUN,KNUL,ISTAT)   ! Rotation angle
      CALL STKRDR('INPUTR',2,1,IAV,RBIN,KUN,KNUL,ISTAT) ! Binning factor
      PI     = 3.14159
      A      = OUTANG*PI/180.
      BIN    = RBIN
      ERROR  = 0
      IF (METHOD(1:1) .EQ. 'p')  WRITE(METHOD,111) 'PAIR'
      IF (METHOD(1:1) .EQ. 'a')  WRITE(METHOD,111) 'ANGLE'
111   FORMAT(A5)

      N = 0
      DO 100 ROW = NROW1, 1, -1
        CALL TBSGET(TID, ROW, SELECT, ISTAT)
        IF (SELECT) THEN
            N = N + 1
            CALL TBERDD(TID, ROW, ICOLO,  IP(N), KNUL, ISTAT)
            CALL TBERDD(TID, ROW, ICOLX,  IX(N), KNUL, ISTAT)
            CALL TBERDD(TID, ROW, ICOLN,  IY(N), KNUL, ISTAT)
            CALL TBERDD(TID, ROW, ICOLW,  IL(N), KNUL, ISTAT)
        ENDIF
100   CONTINUE
      M = IP(1)
      DO 110 ROW = 1, N
         IP(ROW) = IP(ROW) - M
110   CONTINUE

C Check maximum size of arrays

      IF (MAXDIM .LE. (N+10)) THEN
         CALL STTPUT(
     +   'Buffer size MAXDIM too small in module wlcoptim',ISTAT)
         GOTO 1900
      ENDIF

C Check organisation of data if method PAIR

      IF (METHOD(1:1) .EQ. 'P') THEN

         IF (IL(1) .NE. IL(2))     ERROR = 1
         IF (IL(3) .NE. IL(4))     ERROR = 1
         IF (IP(2) .NE. (IP(1)+1)) ERROR = 1
         IF (IP(4) .NE. (IP(3)+1)) ERROR = 1
         IF (N   .NE. 4)           ERROR = 1

         IF (ERROR .EQ. 1) THEN

           IF (IL(1).EQ.IL(3).AND.IL(2).EQ.IL(4).AND.IP(3).EQ.
     C       (IP(1)+1).AND.IP(4).EQ.(IP(2)+1).AND.N.EQ.4) THEN

              ERROR  = 0

              BUFFER = IL(2)
              IL(2)  = IL(3)
              IL(3)  = BUFFER

              BUFFER = IX(2)
              IX(2)  = IX(3)
              IX(3)  = BUFFER

              BUFFER = IY(2)
              IY(2)  = IY(3)
              IY(3)  = BUFFER

              BUFFER = IP(2)
              IP(2)  = IP(3)
              IP(3)  = NINT(BUFFER)

           ELSE
         CALL STTPUT(
     +   'Error:     Identifications are wrongly organized',ISTAT)
         CALL STTPUT(
     +   'Remember : First line:  left-hand side, then right-side',
     +   ISTAT)
         CALL STTPUT(
     +   '           Second line: left-hand side, then right-side',
     +   ISTAT)
         GOTO 1900

           ENDIF

         ENDIF

       ENDIF

C Check that identifications are in different orders for method ANGLE

       IF (METHOD(1:1) .EQ. 'A') THEN

          SPL = 0.D0
          DO 200 ROW = 1,N
             SPL = SPL + IP(ROW)
200       CONTINUE
C234567890123456789012345678901234567890123456789012345678901234567890123456
          IF (SPL .LT. 0.5) THEN
             CALL STTPUT(
     + 'Error: Identifications cannot be all in the same order',ISTAT)
             GOTO 1900
          ENDIF

       ENDIF

C Computes sums

        SX    = 0.D0
        SY    = 0.D0
        SX2   = 0.D0
        SXY   = 0.D0
        SY2   = 0.D0

        SYL   = 0.D0
        SYPL  = 0.D0
        SPL   = 0.D0
        SL    = 0.D0
        SL2   = 0.D0
        SPL2  = 0.D0
        SP2L2 = 0.D0
        SXL   = 0.D0
        SXPL  = 0.D0

        DO 10 LOOP = 1 , N

           SX    = SX    + IX(LOOP)
           SY    = SY    + IY(LOOP)
           SX2   = SX2   + IX(LOOP)*IX(LOOP)
           SY2   = SY2   + IY(LOOP)*IY(LOOP)
           SXY   = SXY   + IX(LOOP)*IY(LOOP)
           SYL   = SYL   + IY(LOOP)*IL(LOOP)
           SYPL  = SYPL  + IY(LOOP)*IP(LOOP)*IL(LOOP)
           SPL   = SPL   + IP(LOOP)*IL(LOOP)
           SL    = SL    + IL(LOOP)
           SL2   = SL2   + IL(LOOP)*IL(LOOP)
           SPL2  = SPL2  + IP(LOOP)*IL(LOOP)**2
           SP2L2 = SP2L2 + IP(LOOP)**2*IL(LOOP)**2
           SXL   = SXL   + IX(LOOP)*IL(LOOP)
           SXPL  = SXPL  + IX(LOOP)*IP(LOOP)*IL(LOOP)

10      CONTINUE

      SXML = M*SXL + SXPL
      SYML = M*SYL + SYPL
      SML  = M*SL  + SPL
      SML2 = M*M*SL2 + 2*M*SPL2 + SP2L2

C  Computes a linear fit and display results

      SXR   = COS(A)*SX + BIN*SIN(A)*SY
      SXR2  = SX2*COS(A)**2+2*SXY*COS(A)*BIN*SIN(A)+
     1        SY2*SIN(A)**2*BIN**2
      SXRML = SXML*COS(A)+SYML*SIN(A)*BIN
      DET   = N*SXR2-SXR*SXR
      COEFA = (SML*SXR2-SXR*SXRML)/DET
      COEFB = (N*SXRML-SXR*SML)/DET
      RMS   = SML2-2*COEFA*COEFB*SXR-COEFB*COEFB*SXR2-N*COEFA*COEFA
      RMS   = SQRT(RMS/N)

C     Display header

      WRITE(TEXT,510)
      CALL STTPUT(TEXT,ISTAT)
      WRITE(TEXT,570)
      CALL STTPUT(TEXT,ISTAT)
      WRITE(TEXT,510)
      CALL STTPUT(TEXT,ISTAT)
      WRITE(TEXT,520)
      CALL STTPUT(TEXT,ISTAT)
      WRITE(TEXT,500) COEFA,COEFB,RMS
      CALL STTPUT(TEXT,ISTAT)
      WRITE(TEXT,510)
      CALL STTPUT(TEXT,ISTAT)
      WRITE(TEXT,530)
      CALL STTPUT(TEXT,ISTAT)
      WRITE(TEXT,510)
      CALL STTPUT(TEXT,ISTAT)

500   FORMAT (1X,'* A = ',F12.5,' B= ',F12.5,' RMS = ',F12.3,7X,'*')
510   FORMAT (1X,61('*'))
520   FORMAT (1X,'*',10X,'Linear fit of : m * lambda = b * xr + a',
     C        10X,'*')
530   FORMAT (1X,'*    Parameter    *    Value    *    Optimum   *',
     C        '    RMS     *')
540   FORMAT (1X,'*',A16,' *',F12.5,' *',F13.5,' *',F11.2,' *')
550   FORMAT (1X,'*',' Ident ',I6,4X,'*',F12.5,' *',F13.5,' *',
     C        F11.2,' *')
560   FORMAT (1X,'*',A16,' *   ',I6,'    *   ',I7,'    *',F11.2,' *')
570   FORMAT (1X,'*',13X,'Wavelength Calibration Diagnosis',14X,'*')
C Least square solutions for the methods PAIR and ANGLE

C   Optimal Rotation Angle

      DET = ABS(N*SYML-SML*SY)

      IF (DET .GT. 0.01) THEN

      TGA1 = (N*(SX2*SYML-SXML*SXY)-SX**2*SYML+SML*(SX*SXY-SX2
     1   *SY)+SX*SXML*SY)/(BIN*(N*(SXML*SY2-SXY*SYML)+SX*SY*SYML+SML*(SX
     2   Y*SY-SX*SY2)-SXML*SY**2))

      TGA2 = -(N*SXML-SML*SX)/(BIN*(N*SYML-SML*SY))

      AOPT = ATAN(TGA1)

      ELSE

      AOPT = -999
   
      ENDIF

C   Optimal Binning Factor

      IF (ABS(A) .GT. 0.02) THEN

      B1 = COS(A)*(N*(SX2*SYML-SXML*SXY)-SX**2*SYML+SML*(SX*SXY-SX2*SY
     1   )+SX*SXML*SY)/(SIN(A)*(N*(SXML*SY2-SXY*SYML)+SX*SY*SYML+SML*(SX
     2   Y*SY-SX*SY2)-SXML*SY**2))

      ELSE

      B1 = -999

      ENDIF

C   Optimal Absolute Order Number for the first identified line

      MDP = -((SIN(A)**2*BIN**2*N*SYL-SIN(A)**2*BIN**2*SL*SY+COS(
     1   A)*SIN(A)*BIN*N*SXL-COS(A)*SIN(A)*BIN*SL*SX)*SYPL+(-SIN(A)**2*B
     2   IN**2*SPL*SY+COS(A)*SIN(A)*BIN*N*SXPL-COS(A)*SIN(A)*BIN*SPL*SX)
     3   *SYL+(SIN(A)**2*BIN**2*SL*SPL-SIN(A)**2*BIN**2*N*SPL2)*SY2+SIN(
     4   A)**2*BIN**2*SPL2*SY**2+(-COS(A)*SIN(A)*BIN*SL*SXPL-COS(A)*SIN(
     5   A)*BIN*SPL*SXL+2*COS(A)*SIN(A)*BIN*SPL2*SX)*SY+(2*COS(A)*SIN(A)
     6   *BIN*SL*SPL-2*COS(A)*SIN(A)*BIN*N*SPL2)*SXY+(COS(A)**2*N*SXL-CO
     7   S(A)**2*SL*SX)*SXPL-COS(A)**2*SPL*SX*SXL+(COS(A)**2*SL*SPL-COS(
     8   A)**2*N*SPL2)*SX2+COS(A)**2*SPL2*SX**2)/(SIN(A)**2*BIN**2*N*SYL
     9   **2+(-2*SIN(A)**2*BIN**2*SL*SY+2*COS(A)*SIN(A)*BIN*N*SXL-2*COS(
     :   A)*SIN(A)*BIN*SL*SX)*SYL+(SIN(A)**2*BIN**2*SL**2-SIN(A)**2*BIN*
     ;   *2*N*SL2)*SY2+SIN(A)**2*BIN**2*SL2*SY**2+(2*COS(A)*SIN(A)*BIN*S
     <   L2*SX-2*COS(A)*SIN(A)*BIN*SL*SXL)*SY+(2*COS(A)*SIN(A)*BIN*SL**2
     =   -2*COS(A)*SIN(A)*BIN*N*SL2)*SXY+COS(A)**2*N*SXL**2-2*COS(A)**2*
     >   SL*SX*SXL+(COS(A)**2*SL**2-COS(A)**2*N*SL2)*SX2+COS(A)**2*SL2*S
     ?   X**2)

      M2OPT = NINT(MDP)

C Compute new rms if a value is changed.

      MOLD = M         ! Keep the old value in memory
      M = M2OPT        ! Replace current value by optimum
      MODE = 'M'       ! Set the mode
      GOTO 1000        ! Compute rms 
1001  CONTINUE         ! Return point
      M = MOLD         ! Restore old value
      MRMS = RMSM      ! Store rms
      MINRMS(1) = RMSM

      WRITE(TEXT,560) 'Order number',M,M2OPT,MRMS
      CALL STTPUT(TEXT,ISTAT)
      CALL STKWRI('OUTPUTI',M2OPT,1,1,KUN,ISTAT)

      AOLD = A
      A = AOPT
      MODE = 'A'
      GOTO 1000
1002  CONTINUE         ! Return point
      A = AOLD
      ARMS = RMSM
      MINRMS(2) = RMSM

      WRITE(TEXT,540) 'Rotation angle',(A*180./PI),(AOPT*180/PI),ARMS
      CALL STTPUT(TEXT,ISTAT)
      OUTANG = AOPT*180./PI
      CALL STKWRR('OUTPUTR',OUTANG,1,1,KUN,ISTAT)

      BOLD = BIN
      BIN = B1
      MODE = 'B'
      GOTO 1000
1003  CONTINUE         ! Return point
      BIN = BOLD
      BRMS = RMSM
      MINRMS(3) = RMSM

      IF (B1 .GT. 0) THEN
        WRITE(TEXT,540) 'Binning factor',BIN,SQRT(B1),BRMS
        CALL STTPUT(TEXT,ISTAT)
      ELSE
        WRITE(TEXT,540) 'Binning factor',BIN,BIN,0.
        CALL STTPUT(TEXT,ISTAT)
      ENDIF

C
C Least square solutions for the method PAIR

      IF (METHOD(1:1) .EQ. 'P') THEN

      X1  = IX(1)
      XP1 = IX(2)
      X2  = IX(3)
      XP2 = IX(4)
      P   = IP(3) - IP(1)
      L1  = IL(1)
      L2  = IL(3)

      L1OPT = (((2*L2*M+L2)*P+2*L2*M**2+L2*M)*XP2**2+((L2*P+2*L2*M+2*L2)
     1   *XP1+((-4*L2*M-2*L2)*P-4*L2*M**2-4*L2*M-L2)*X2+(-L2*P-L2)*X1)
     2   *XP2+(2*L2*M*P+2*L2*M**2+L2*M)*XP1**2+((L2*P-L2)*X2+
     3   ((-4*L2*M-2*L2)*P-4*L2*M**2-4*L2*M-L2)*X1)*XP1+((2*L2*M+L2)*P+
     4   2*L2*M**2+3*L2*M+L2)*X2**2+(-L2*P-2*L2*M)*X1*X2+((2*L2*M+2*L2)
     5   *P+2*L2*M**2+3*L2*M+L2)*X1**2)/((2*M**2+2*M+2)*XP2**2+(2*M*XP1
     6   +(-4*M**2-4*M-2)*X2+(-2*M-2)*X1)*XP2+2*M**2*XP1**2+(2*M*X2+(-4
     7   *M**2-4*M)*X1)*XP1+(2*M**2+2*M+2)*X2**2+(-2*M-2)*X1*X2+(2*M**2
     8   +4*M+2)*X1**2)

         L1OLD = L1
         L1 = L1OPT

      RMSL1 = -4*((2*L2*P+2*(L2+L1)*M+L2+L1)*(XP2**2+XP1**2+X2**2+X1**2)
     1 -(XP2+XP1+X2+X1)*(L2*(P+M+1)*XP2+L1*(M+1)*XP1+L2*(P+M)*X2+L1*M*X1
     2   ))**2/(4*(XP2**2+XP1**2+X2**2+X1**2)-(XP2+XP1+X2+X1)**2)**2-2*(
     3   XP2+XP1+X2+X1)*(4*(L2*(P+M+1)*XP2+L1*(M+1)*XP1+L2*(P+M)*X2+L1*M
     4   *X1)-(2*L2*P+2*(L2+L1)*M+L2+L1)*(XP2+XP1+X2+X1))*((2*L2*P+2*(L2
     5   +L1)*M+L2+L1)*(XP2**2+XP1**2+X2**2+X1**2)-(XP2+XP1+X2+X1)*(L2*(
     6   P+M+1)*XP2+L1*(M+1)*XP1+L2*(P+M)*X2+L1*M*X1))/(4*(XP2**2+XP1**2
     7   +X2**2+X1**2)-(XP2+XP1+X2+X1)**2)**2-(4*(L2*(P+M+1)*XP2+L1*(M+1
     8   )*XP1+L2*(P+M)*X2+L1*M*X1)-(2*L2*P+2*(L2+L1)*M+L2+L1)*(XP2+XP1+
     9   X2+X1))**2*(XP2**2+XP1**2+X2**2+X1**2)/(4*(XP2**2+XP1**2+X2**2+
     :   X1**2)-(XP2+XP1+X2+X1)**2)**2+L2**2*(P+M+1)**2+L2**2*(P+M)**2+L
     ;   1**2*(M+1)**2+L1**2*M**2

         L1 = L1OLD

      L2OPT = (((2*L1*M+L1)*P+2*L1*M**2+L1*M)*XP2**2+((L1*P+2*L1*M+2
     1   *L1)*XP1+((-4*L1*M-2*L1)*P-4*L1*M**2-4*L1*M-L1)*X2+(-L1*P-L1)*X
     2   1)*XP2+(2*L1*M*P+2*L1*M**2+L1*M)*XP1**2+((L1*P-L1)*X2+((-4*L1*M
     3   -2*L1)*P-4*L1*M**2-4*L1*M-L1)*X1)*XP1+((2*L1*M+L1)*P+2*L1*M**2+
     4   3*L1*M+L1)*X2**2+(-L1*P-2*L1*M)*X1*X2+((2*L1*M+2*L1)*P+2*L1*M**
     5   2+3*L1*M+L1)*X1**2)/((2*P**2+4*M*P+2*M**2)*XP2**2+((2*P+2*M)*XP
     6   1+(-4*P**2+(-8*M-4)*P-4*M**2-4*M)*X2+(2*P+2*M)*X1)*XP2+(2*P**2+
     7   (4*M+2)*P+2*M**2+2*M+2)*XP1**2+((-2*P-2*M-2)*X2+(-4*P**2+(-8*M-
     8   4)*P-4*M**2-4*M-2)*X1)*XP1+(2*P**2+(4*M+4)*P+2*M**2+4*M+2)*X2**
     9   2+(-2*P-2*M-2)*X1*X2+(2*P**2+(4*M+2)*P+2*M**2+2*M+2)*X1**2)

         L2OLD = L2
         L2 = L2OPT

      RMSL2 = -4*((2*L2*P+2*(L2+L1)*M+L2+L1)*(XP2**2+XP1**2+X2**2+X1**2)
     1 -(XP2+XP1+X2+X1)*(L2*(P+M+1)*XP2+L1*(M+1)*XP1+L2*(P+M)*X2+L1*M*X1
     2   ))**2/(4*(XP2**2+XP1**2+X2**2+X1**2)-(XP2+XP1+X2+X1)**2)**2-2*(
     3   XP2+XP1+X2+X1)*(4*(L2*(P+M+1)*XP2+L1*(M+1)*XP1+L2*(P+M)*X2+L1*M
     4   *X1)-(2*L2*P+2*(L2+L1)*M+L2+L1)*(XP2+XP1+X2+X1))*((2*L2*P+2*(L2
     5   +L1)*M+L2+L1)*(XP2**2+XP1**2+X2**2+X1**2)-(XP2+XP1+X2+X1)*(L2*(
     6   P+M+1)*XP2+L1*(M+1)*XP1+L2*(P+M)*X2+L1*M*X1))/(4*(XP2**2+XP1**2
     7   +X2**2+X1**2)-(XP2+XP1+X2+X1)**2)**2-(4*(L2*(P+M+1)*XP2+L1*(M+1
     8   )*XP1+L2*(P+M)*X2+L1*M*X1)-(2*L2*P+2*(L2+L1)*M+L2+L1)*(XP2+XP1+
     9   X2+X1))**2*(XP2**2+XP1**2+X2**2+X1**2)/(4*(XP2**2+XP1**2+X2**2+
     :   X1**2)-(XP2+XP1+X2+X1)**2)**2+L2**2*(P+M+1)**2+L2**2*(P+M)**2+L
     ;   1**2*(M+1)**2+L1**2*M**2


         RMSL1 = SQRT(RMSL1/4.)
         RMSL2 = SQRT(RMSL2/4.)

        WRITE(TEXT,540) 'Ident. 1 and 2',L1OLD,L1OPT,RMSL1
        CALL STTPUT(TEXT,ISTAT)
        WRITE(TEXT,540) 'Ident. 3 and 4',L2OLD,L2OPT,RMSL2
        CALL STTPUT(TEXT,ISTAT)

        MINRMS(4) = RMSL1
        MINRMS(5) = RMSL1
        MINRMS(6) = RMSL2
        MINRMS(7) = RMSL2
        LKOPT(1)  = L1OPT
        LKOPT(2)  = L1OPT
        LKOPT(3)  = L2OPT
        LKOPT(4)  = L2OPT

      ENDIF


C  Optimal Wavelength for each of the computed wavelength IL(LOOP)

      IF (METHOD(1:1) .EQ. 'A') THEN


      DO 50 K = 1,N

        PK  = IP(K)
        LK  = IL(K)
        XK  = IX(K)
        YK  = IY(K)

        SYLK   = 0.D0
        SYPLK  = 0.D0
        SPLK   = 0.D0
        SLK    = 0.D0
        SL2K   = 0.D0
        SPL2K  = 0.D0
        SP2L2K = 0.D0
        SXLK   = 0.D0
        SXPLK  = 0.D0

        DO 30 LOOP = 1 , N

         IF (LOOP.NE.K) THEN
           SYLK   = SYLK   + IY(LOOP)*IL(LOOP)
           SYPLK  = SYPLK  + IY(LOOP)*IP(LOOP)*IL(LOOP)
           SPLK   = SPLK   + IP(LOOP)*IL(LOOP)
           SLK    = SLK    + IL(LOOP)
           SL2K   = SL2K   + IL(LOOP)*IL(LOOP)
           SPL2K  = SPL2K  + IP(LOOP)*IL(LOOP)**2
           SP2L2K = SP2L2K + IP(LOOP)**2*IL(LOOP)**2
           SXLK   = SXLK   + IX(LOOP)*IL(LOOP)
           SXPLK  = SXPLK  + IX(LOOP)*IP(LOOP)*IL(LOOP)
         ENDIF

30      CONTINUE

C The formula for LKOPT is splitted in two parts (too big for the poor
C  compiler !!)

        SBIN = (
     <   (SIN(A)**2*BIN**2*N*PK+SIN(A)**2*BIN**2*M*N)*YK**2+((2*COS(A)*S
     =   IN(A)*BIN*N*PK+2*COS(A)*SIN(A)*BIN*M*N)*XK+(-2*SIN(A)**2*BIN**2
     >   *PK-2*SIN(A)**2*BIN**2*M)*SY+(-2*COS(A)*SIN(A)*BIN*PK-2*COS(A)*
     ?   SIN(A)*BIN*M)*SX)*YK+(COS(A)**2*N*PK+COS(A)**2*M*N)*XK**2+((-2*
     @   COS(A)*SIN(A)*BIN*PK-2*COS(A)*SIN(A)*BIN*M)*SY+(-2*COS(A)**2*PK
     1   -2*COS(A)**2*M)*SX)*XK+((SIN(A)**2*BIN**2-SIN(A)**2*BIN**2*N)*P
     2   K-SIN(A)**2*BIN**2*M*N+SIN(A)**2*BIN**2*M)*SY2+(SIN(A)**2*BIN**
     3   2*PK+SIN(A)**2*BIN**2*M)*SY**2+(2*COS(A)*SIN(A)*BIN*PK+2*COS(A)
     4   *SIN(A)*BIN*M)*SX*SY+((2*COS(A)*SIN(A)*BIN-2*COS(A)*SIN(A)*BIN*
     5   N)*PK-2*COS(A)*SIN(A)*BIN*M*N+2*COS(A)*SIN(A)*BIN*M)*SXY+((COS(
     6   A)**2-COS(A)**2*N)*PK-COS(A)**2*M*N+COS(A)**2*M)*SX2+(COS(A)**2
     7   *PK+COS(A)**2*M)*SX**2)

       LKOPT(K) = -((SIN(A)**2*BIN**2*N*SYPLK+SIN(A)**2*BIN**2*M*N*SY
     1   LK+(-SIN(A)**2*BIN**2*SPLK-SIN(A)**2*BIN**2*M*SLK)*SY+COS(A)*SI
     2   N(A)*BIN*N*SXPLK+COS(A)*SIN(A)*BIN*M*N*SXLK+(-COS(A)*SIN(A)*BIN
     3   *SPLK-COS(A)*SIN(A)*BIN*M*SLK)*SX)*YK+(COS(A)*SIN(A)*BIN*N*SYPL
     4   K+COS(A)*SIN(A)*BIN*M*N*SYLK+(-COS(A)*SIN(A)*BIN*SPLK-COS(A)*SI
     5   N(A)*BIN*M*SLK)*SY+COS(A)**2*N*SXPLK+COS(A)**2*M*N*SXLK+(-COS(A
     6   )**2*SPLK-COS(A)**2*M*SLK)*SX)*XK+(-SIN(A)**2*BIN**2*SY-COS(A)*
     7   SIN(A)*BIN*SX)*SYPLK+(-SIN(A)**2*BIN**2*M*SY-COS(A)*SIN(A)*BIN*
     8   M*SX)*SYLK+(SIN(A)**2*BIN**2*SPLK+SIN(A)**2*BIN**2*M*SLK)*SY2+(
     9   -COS(A)*SIN(A)*BIN*SXPLK-COS(A)*SIN(A)*BIN*M*SXLK)*SY+(2*COS(A)
     :   *SIN(A)*BIN*SPLK+2*COS(A)*SIN(A)*BIN*M*SLK)*SXY-COS(A)**2*SX*SX
     ;   PLK-COS(A)**2*M*SX*SXLK+(COS(A)**2*SPLK+COS(A)**2*M*SLK)*SX2)/
     +   SBIN

         LKO = LKOPT(K)

      RMSLK = -N*((SPLK+M*SLK+LKO*(PK+M))*(SIN(A)**2*BIN**2*SY2+2*COS(A)
     1   *SIN(A)*BIN*SXY+COS(A)**2*SX2)-(SIN(A)*BIN*SY+COS(A)*SX)*(SIN(A
     2   )*BIN*(M*(LKO*YK+SYLK)+LKO*PK*YK+SYPLK)+COS(A)*(M*(LKO*XK+SXLK)
     C   +LKO
     3   *PK*XK+SXPLK)))**2/(N*(SIN(A)**2*BIN**2*SY2+2*COS(A)*SIN(A)*BIN
     4   *SXY+COS(A)**2*SX2)-(SIN(A)*BIN*SY+COS(A)*SX)**2)**2-2*(SIN(A)*
     5 BIN*SY+COS(A)*SX)*(N*(SIN(A)*BIN*(M*(LKO*YK+SYLK)+LKO*PK*YK+SYPLK
     6)+COS(A)*(M*(LKO*XK+SXLK)+LKO*PK*XK+SXPLK))-(SPLK+M*SLK+LKO*(PK+M)
     7  )*(SIN(A)*BIN*SY+COS(A)*SX))*((SPLK+M*SLK+LKO*(PK+M))*(SIN(A)**2
     8   *BIN**2*SY2+2*COS(A)*SIN(A)*BIN*SXY+COS(A)**2*SX2)-(SIN(A)*BIN*
     9 SY+COS(A)*SX)*(SIN(A)*BIN*(M*(LKO*YK+SYLK)+LKO*PK*YK+SYPLK)+COS(A
     : )*(M*(LKO*XK+SXLK)+LKO*PK*XK+SXPLK)))/(N*(SIN(A)**2*BIN**2*SY2+2*
     ;   COS(A)*SIN(A)*BIN*SXY+COS(A)**2*SX2)-(SIN(A)*BIN*SY+COS(A)*SX)*
     <   *2)**2-(SIN(A)**2*BIN**2*SY2+2*COS(A)*SIN(A)*BIN*SXY+COS(A)**2*
     + SX2)*(N*(SIN(A)*BIN*(M*(LKO*YK+SYLK)+LKO*PK*YK+SYPLK)+COS(A)*(M*(
     >LKO*XK+SXLK)+LKO*PK*XK+SXPLK))-(SPLK+M*SLK+LKO*(PK+M))*(SIN(A)*BIN
     ?   *SY+COS(A)*SX))**2/(N*(SIN(A)**2*BIN**2*SY2+2*COS(A)*SIN(A)*BIN
     @   *SXY+COS(A)**2*SX2)-(SIN(A)*BIN*SY+COS(A)*SX)**2)**2+2*M*SPL2K+
     1   SP2L2K+M**2*SL2K+LKO**2*(PK+M)**2


         RMSLK = SQRT(RMSLK/N)
         MINRMS(3+K) = RMSLK
         WRITE(TEXT,550) K,LK,LKOPT(K),RMSLK
         CALL STTPUT(TEXT,ISTAT)

50      CONTINUE

        ENDIF

        WRITE(TEXT,510)
        CALL STTPUT(TEXT,ISTAT)

        GOTO 2000

C **********************************************************************

1000    CONTINUE
C       Computes RMS(A,M,BIN). There are too many parameters to create
C       a subroutine. Let's use the good old GOTO.

       RMSM = -N*((SPL+M*SL)*(SIN(A)**2*BIN**2*SY2+2*COS(A)*SIN(A)*BIN*S
     1   XY+COS(A)**2*SX2)-(SIN(A)*BIN*SY+COS(A)*SX)*(SIN(A)*BIN*(SYPL+M
     2   *SYL)+COS(A)*(SXPL+M*SXL)))**2/(N*(SIN(A)**2*BIN**2*SY2+2*COS(A
     3   )*SIN(A)*BIN*SXY+COS(A)**2*SX2)-(SIN(A)*BIN*SY+COS(A)*SX)**2)**
     4   2-2*(SIN(A)*BIN*SY+COS(A)*SX)*(N*(SIN(A)*BIN*(SYPL+M*SYL)+COS(A
     5   )*(SXPL+M*SXL))-(SPL+M*SL)*(SIN(A)*BIN*SY+COS(A)*SX))*((SPL+M*S
     6   L)*(SIN(A)**2*BIN**2*SY2+2*COS(A)*SIN(A)*BIN*SXY+COS(A)**2*SX2)
     7   -(SIN(A)*BIN*SY+COS(A)*SX)*(SIN(A)*BIN*(SYPL+M*SYL)+COS(A)*(SXP
     8   L+M*SXL)))/(N*(SIN(A)**2*BIN**2*SY2+2*COS(A)*SIN(A)*BIN*SXY+COS
     9   (A)**2*SX2)-(SIN(A)*BIN*SY+COS(A)*SX)**2)**2-(SIN(A)**2*BIN**2*
     :   SY2+2*COS(A)*SIN(A)*BIN*SXY+COS(A)**2*SX2)*(N*(SIN(A)*BIN*(SYPL
     ;   +M*SYL)+COS(A)*(SXPL+M*SXL))-(SPL+M*SL)*(SIN(A)*BIN*SY+COS(A)*S
     <   X))**2/(N*(SIN(A)**2*BIN**2*SY2+2*COS(A)*SIN(A)*BIN*SXY+COS(A)*
     =   *2*SX2)-(SIN(A)*BIN*SY+COS(A)*SX)**2)**2+2*M*SPL2+SP2L2+M**2*SL
     >   2

         RMSM = SQRT(RMSM/N)

        IF (MODE(1:1) .EQ. 'M') GOTO 1001  ! Absolute order number
        IF (MODE(1:1) .EQ. 'A') GOTO 1002  ! Angle
        IF (MODE(1:1) .EQ. 'B') GOTO 1003  ! Bin

C **********************************************************************

1900    CONTINUE
        ERROR = 1

2000    CONTINUE
        CALL STKWRI('OUTPUTI',ERROR,1,1,KUN,ISTAT)
        CALL STSEPI
        STOP
        END





