← Back to branch summary

~ubuntu-branches/debian/jessie/eso-midas/jessie

« back to all changes in this revision

Viewing changes to stdred/echelle/src/necwlcopt.for

  • Committer: Package Import Robot
  • Author(s): Ole Streicher
  • Date: 2014-04-22 14:44:58 UTC
  • Revision ID: package-import@ubuntu.com-20140422144458-okiwi1assxkkiz39
Tags: upstream-13.09pl1.2+dfsg
Import upstream version 13.09pl1.2+dfsg

Show diffs side-by-side

added added

removed removed

Lines of Context:
stdred/echelle/src/necwlcopt.for
 
1
C @(#)necwlcopt.for     19.1 (ES0-DMD) 02/25/03 14:20:26
 
2
C===========================================================================
 
3
C Copyright (C) 1995 European Southern Observatory (ESO)
 
4
C
 
5
C This program is free software; you can redistribute it and/or 
 
6
C modify it under the terms of the GNU General Public License as 
 
7
C published by the Free Software Foundation; either version 2 of 
 
8
C the License, or (at your option) any later version.
 
9
C
 
10
C This program is distributed in the hope that it will be useful,
 
11
C but WITHOUT ANY WARRANTY; without even the implied warranty of
 
12
C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
13
C GNU General Public License for more details.
 
14
C
 
15
C You should have received a copy of the GNU General Public 
 
16
C License along with this program; if not, write to the Free 
 
17
C Software Foundation, Inc., 675 Massachusetss Ave, Cambridge, 
 
18
C MA 02139, USA.
 
19
C
 
20
C Corresponding concerning ESO-MIDAS should be addressed as follows:
 
21
C       Internet e-mail: midas@eso.org
 
22
C       Postal address: European Southern Observatory
 
23
C                       Data Management Division 
 
24
C                       Karl-Schwarzschild-Strasse 2
 
25
C                       D 85748 Garching bei Muenchen 
 
26
C                       GERMANY
 
27
C===========================================================================
 
28
C
 
29
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
30
C.COPYRIGHT   (C) 1992 European Southern Observatory
 
31
C.IDENT       wlcoptim.for
 
32
C.AUTHOR      Pascal Ballester,  ESO - Garching
 
33
C.KEYWORDS    Spectroscopy, Echelle,
 
34
C.PURPOSE     Computes optimal parameters for the 
 
35
C             wavelength calibration
 
36
C.VERSION     1.0    Creation    1992-MAY-26
 
37
C-------------------------------------------------------
 
38
C
 
39
C
 
40
C  Program : Wavelength Calibration Least Squares Solutions
 
41
C
 
42
C  This program has been partly generated using Macsyma to 
 
43
C  solve the following least-square problem:
 
44
C
 
45
C     Xr(i) = cos(a) * X(i) + bin * sin(a) * Y(i)    (1)
 
46
C
 
47
C  takes into account a binning factor of the ccd and a rotation
 
48
C  of the detector.
 
49
C
 
50
C     F(Xr(i)) = B * Xr(i) + A = (m+p(i))*lambda(i)  (2)
 
51
C
 
52
C  is the echelle relation applied in the rotated space and approximated
 
53
C  by a linear relation. m is the absolute order number of a line
 
54
C  and p(i) is relative order number of the other lines.
 
55
C
 
56
C  Macsyma has been used to find different optima, like:
 
57
C    - The absolute order number of the first identification : m
 
58
C    - The binning factor    : bin
 
59
C    - In method ANGLE, the angle of rotation : a
 
60
C    - In method PAIR, the wavelength of the two pairs : lambda(1), lambda(2)
 
61
C
 
62
C  The different least-square solutions are used to compute the optimal
 
63
C  value of a parameter, assuming that all other parameters are correctly
 
64
C  provided. This is a cross-check of the entered values. By computing 
 
65
C  the new rms after replacement of a parameter by its optimal value,
 
66
C  it is possible to provide a diagnosis of what action would benefit
 
67
C  the most to solve the problem.
 
68
C
 
69
C  INPUT:
 
70
C   IN_A/C/1/60   Name of line table
 
71
C   IN_B/C/1/60   Method (PAIR or ANGLE)
 
72
C   INPUTR/R/1/2  Angle of rotation, Binning ratio.
 
73
C
 
74
C  OUTPUT:
 
75
C   OUTPUTI/I/1/1 Order number of the first line
 
76
C   OUTPUTR/R/1/2 Optimal rotation angle, Optimal binning factor
 
77
C
 
78
        IMPLICIT NONE
 
79
 
 
80
        INTEGER  MAXDIM
 
81
        PARAMETER (MAXDIM=1500)
 
82
        DOUBLE PRECISION    IX(MAXDIM), IY(MAXDIM)
 
83
        DOUBLE PRECISION    IP(MAXDIM), IL(MAXDIM)
 
84
        DOUBLE PRECISION    SX, SY, SX2, SY2, SXY
 
85
C        DOUBLE PRECISION    SML2K, SXMLK, SYMLK, SMLK
 
86
        DOUBLE PRECISION    SP2L2K, BUFFER
 
87
        DOUBLE PRECISION    SXLK, SXPLK, SYLK, SYPLK, SPLK, SLK
 
88
        DOUBLE PRECISION    SL2K, SPL2K, PK, LK, XK, YK, RMSLK
 
89
        DOUBLE PRECISION    SML2, SP2L2, SXML, SYML, SML
 
90
        DOUBLE PRECISION    SXL, SXPL, SYL, SYPL, SPL, SL
 
91
        DOUBLE PRECISION    SL2, SPL2, SBIN
 
92
        INTEGER             M, N, K, M2OPT, LOOP, ERROR, P
 
93
        DOUBLE PRECISION    BIN, PI, MDP, MINRMS(MAXDIM)
 
94
        DOUBLE PRECISION    A, B1
 
95
        DOUBLE PRECISION    TGA1, TGA2, AOPT
 
96
        REAL                OUTANG, RBIN
 
97
        DOUBLE PRECISION    L1,L2,X1,X2,XP1,XP2,L1OPT,L2OPT
 
98
        DOUBLE PRECISION    SXR, SXR2, SXRML, DET, COEFA, COEFB, RMS
 
99
        DOUBLE PRECISION    L1OLD, L2OLD, RMSL1, RMSL2
 
100
        CHARACTER           MODE*5, METHOD*10, TEXT*80
 
101
        INTEGER             MOLD, ROW
 
102
        DOUBLE PRECISION    RMSM, MRMS, ARMS, BRMS, AOLD, BOLD
 
103
        DOUBLE PRECISION    LKOPT(MAXDIM), LKO
 
104
 
 
105
        CHARACTER           TABLE*60
 
106
        INTEGER             IAV, KUN, KNUL, ISTAT, ICOLX, ICOLW
 
107
        INTEGER             ICOLO, ICOLN, TID, NCOL1, NROW1, NSC1
 
108
        INTEGER             NACOL, NAROW, MADRID
 
109
        LOGICAL             SELECT
 
110
 
 
111
        INCLUDE 'MID_INCLUDE:ST_DEF.INC'
 
112
        COMMON/VMR/MADRID(1)
 
113
        INCLUDE 'MID_INCLUDE:ST_DAT.INC'
 
114
 
 
115
C
 
116
C Initialisations
 
117
 
 
118
C
 
119
C ... Init input table
 
120
C
 
121
      CALL STSPRO('WLCLS')
 
122
 
 
123
      CALL STKRDC('IN_A',1,1,60,IAV,TABLE,KUN,KNUL,ISTAT)
 
124
      CALL TBTOPN(TABLE,F_U_MODE,TID,ISTAT)
 
125
C
 
126
      CALL TBIGET(TID,NCOL1,NROW1,NSC1,NACOL,NAROW,ISTAT)
 
127
      CALL TBCSER(TID,':X '    ,ICOLX,ISTAT)
 
128
      CALL TBCSER(TID,':IDENT' ,ICOLW,ISTAT)
 
129
      CALL TBCSER(TID,':ORDER' ,ICOLO,ISTAT)
 
130
      CALL TBCSER(TID,':YNEW'  ,ICOLN,ISTAT)
 
131
 
 
132
      CALL STKRDC('IN_B',1,1,10,IAV,METHOD,KUN,KNUL,ISTAT)
 
133
      CALL STKRDR('INPUTR',1,1,IAV,OUTANG,KUN,KNUL,ISTAT)   ! Rotation angle
 
134
      CALL STKRDR('INPUTR',2,1,IAV,RBIN,KUN,KNUL,ISTAT) ! Binning factor
 
135
      PI     = 3.14159
 
136
      A      = OUTANG*PI/180.
 
137
      BIN    = RBIN
 
138
      ERROR  = 0
 
139
      IF (METHOD(1:1) .EQ. 'p')  WRITE(METHOD,111) 'PAIR'
 
140
      IF (METHOD(1:1) .EQ. 'a')  WRITE(METHOD,111) 'ANGLE'
 
141
111   FORMAT(A5)
 
142
 
 
143
      N = 0
 
144
      DO 100 ROW = NROW1, 1, -1
 
145
        CALL TBSGET(TID, ROW, SELECT, ISTAT)
 
146
        IF (SELECT) THEN
 
147
            N = N + 1
 
148
            CALL TBERDD(TID, ROW, ICOLO,  IP(N), KNUL, ISTAT)
 
149
            CALL TBERDD(TID, ROW, ICOLX,  IX(N), KNUL, ISTAT)
 
150
            CALL TBERDD(TID, ROW, ICOLN,  IY(N), KNUL, ISTAT)
 
151
            CALL TBERDD(TID, ROW, ICOLW,  IL(N), KNUL, ISTAT)
 
152
        ENDIF
 
153
100   CONTINUE
 
154
      M = IP(1)
 
155
      DO 110 ROW = 1, N
 
156
         IP(ROW) = IP(ROW) - M
 
157
110   CONTINUE
 
158
 
 
159
C Check maximum size of arrays
 
160
 
 
161
      IF (MAXDIM .LE. (N+10)) THEN
 
162
         CALL STTPUT(
 
163
     +   'Buffer size MAXDIM too small in module wlcoptim',ISTAT)
 
164
         GOTO 1900
 
165
      ENDIF
 
166
 
 
167
C Check organisation of data if method PAIR
 
168
 
 
169
      IF (METHOD(1:1) .EQ. 'P') THEN
 
170
 
 
171
         IF (IL(1) .NE. IL(2))     ERROR = 1
 
172
         IF (IL(3) .NE. IL(4))     ERROR = 1
 
173
         IF (IP(2) .NE. (IP(1)+1)) ERROR = 1
 
174
         IF (IP(4) .NE. (IP(3)+1)) ERROR = 1
 
175
         IF (N   .NE. 4)           ERROR = 1
 
176
 
 
177
         IF (ERROR .EQ. 1) THEN
 
178
 
 
179
           IF (IL(1).EQ.IL(3).AND.IL(2).EQ.IL(4).AND.IP(3).EQ.
 
180
     C       (IP(1)+1).AND.IP(4).EQ.(IP(2)+1).AND.N.EQ.4) THEN
 
181
 
 
182
              ERROR  = 0
 
183
 
 
184
              BUFFER = IL(2)
 
185
              IL(2)  = IL(3)
 
186
              IL(3)  = BUFFER
 
187
 
 
188
              BUFFER = IX(2)
 
189
              IX(2)  = IX(3)
 
190
              IX(3)  = BUFFER
 
191
 
 
192
              BUFFER = IY(2)
 
193
              IY(2)  = IY(3)
 
194
              IY(3)  = BUFFER
 
195
 
 
196
              BUFFER = IP(2)
 
197
              IP(2)  = IP(3)
 
198
              IP(3)  = NINT(BUFFER)
 
199
 
 
200
           ELSE
 
201
         CALL STTPUT(
 
202
     +   'Error:     Identifications are wrongly organized',ISTAT)
 
203
         CALL STTPUT(
 
204
     +   'Remember : First line:  left-hand side, then right-side',
 
205
     +   ISTAT)
 
206
         CALL STTPUT(
 
207
     +   '           Second line: left-hand side, then right-side',
 
208
     +   ISTAT)
 
209
         GOTO 1900
 
210
 
 
211
           ENDIF
 
212
 
 
213
         ENDIF
 
214
 
 
215
       ENDIF
 
216
 
 
217
C Check that identifications are in different orders for method ANGLE
 
218
 
 
219
       IF (METHOD(1:1) .EQ. 'A') THEN
 
220
 
 
221
          SPL = 0.D0
 
222
          DO 200 ROW = 1,N
 
223
             SPL = SPL + IP(ROW)
 
224
200       CONTINUE
 
225
C234567890123456789012345678901234567890123456789012345678901234567890123456
 
226
          IF (SPL .LT. 0.5) THEN
 
227
             CALL STTPUT(
 
228
     + 'Error: Identifications cannot be all in the same order',ISTAT)
 
229
             GOTO 1900
 
230
          ENDIF
 
231
 
 
232
       ENDIF
 
233
 
 
234
C Computes sums
 
235
 
 
236
        SX    = 0.D0
 
237
        SY    = 0.D0
 
238
        SX2   = 0.D0
 
239
        SXY   = 0.D0
 
240
        SY2   = 0.D0
 
241
 
 
242
        SYL   = 0.D0
 
243
        SYPL  = 0.D0
 
244
        SPL   = 0.D0
 
245
        SL    = 0.D0
 
246
        SL2   = 0.D0
 
247
        SPL2  = 0.D0
 
248
        SP2L2 = 0.D0
 
249
        SXL   = 0.D0
 
250
        SXPL  = 0.D0
 
251
 
 
252
        DO 10 LOOP = 1 , N
 
253
 
 
254
           SX    = SX    + IX(LOOP)
 
255
           SY    = SY    + IY(LOOP)
 
256
           SX2   = SX2   + IX(LOOP)*IX(LOOP)
 
257
           SY2   = SY2   + IY(LOOP)*IY(LOOP)
 
258
           SXY   = SXY   + IX(LOOP)*IY(LOOP)
 
259
           SYL   = SYL   + IY(LOOP)*IL(LOOP)
 
260
           SYPL  = SYPL  + IY(LOOP)*IP(LOOP)*IL(LOOP)
 
261
           SPL   = SPL   + IP(LOOP)*IL(LOOP)
 
262
           SL    = SL    + IL(LOOP)
 
263
           SL2   = SL2   + IL(LOOP)*IL(LOOP)
 
264
           SPL2  = SPL2  + IP(LOOP)*IL(LOOP)**2
 
265
           SP2L2 = SP2L2 + IP(LOOP)**2*IL(LOOP)**2
 
266
           SXL   = SXL   + IX(LOOP)*IL(LOOP)
 
267
           SXPL  = SXPL  + IX(LOOP)*IP(LOOP)*IL(LOOP)
 
268
 
 
269
10      CONTINUE
 
270
 
 
271
      SXML = M*SXL + SXPL
 
272
      SYML = M*SYL + SYPL
 
273
      SML  = M*SL  + SPL
 
274
      SML2 = M*M*SL2 + 2*M*SPL2 + SP2L2
 
275
 
 
276
C  Computes a linear fit and display results
 
277
 
 
278
      SXR   = COS(A)*SX + BIN*SIN(A)*SY
 
279
      SXR2  = SX2*COS(A)**2+2*SXY*COS(A)*BIN*SIN(A)+
 
280
     1        SY2*SIN(A)**2*BIN**2
 
281
      SXRML = SXML*COS(A)+SYML*SIN(A)*BIN
 
282
      DET   = N*SXR2-SXR*SXR
 
283
      COEFA = (SML*SXR2-SXR*SXRML)/DET
 
284
      COEFB = (N*SXRML-SXR*SML)/DET
 
285
      RMS   = SML2-2*COEFA*COEFB*SXR-COEFB*COEFB*SXR2-N*COEFA*COEFA
 
286
      RMS   = SQRT(RMS/N)
 
287
 
 
288
C     Display header
 
289
 
 
290
      WRITE(TEXT,510)
 
291
      CALL STTPUT(TEXT,ISTAT)
 
292
      WRITE(TEXT,570)
 
293
      CALL STTPUT(TEXT,ISTAT)
 
294
      WRITE(TEXT,510)
 
295
      CALL STTPUT(TEXT,ISTAT)
 
296
      WRITE(TEXT,520)
 
297
      CALL STTPUT(TEXT,ISTAT)
 
298
      WRITE(TEXT,500) COEFA,COEFB,RMS
 
299
      CALL STTPUT(TEXT,ISTAT)
 
300
      WRITE(TEXT,510)
 
301
      CALL STTPUT(TEXT,ISTAT)
 
302
      WRITE(TEXT,530)
 
303
      CALL STTPUT(TEXT,ISTAT)
 
304
      WRITE(TEXT,510)
 
305
      CALL STTPUT(TEXT,ISTAT)
 
306
 
 
307
500   FORMAT (1X,'* A = ',F12.5,' B= ',F12.5,' RMS = ',F12.3,7X,'*')
 
308
510   FORMAT (1X,61('*'))
 
309
520   FORMAT (1X,'*',10X,'Linear fit of : m * lambda = b * xr + a',
 
310
     C        10X,'*')
 
311
530   FORMAT (1X,'*    Parameter    *    Value    *    Optimum   *',
 
312
     C        '    RMS     *')
 
313
540   FORMAT (1X,'*',A16,' *',F12.5,' *',F13.5,' *',F11.2,' *')
 
314
550   FORMAT (1X,'*',' Ident ',I6,4X,'*',F12.5,' *',F13.5,' *',
 
315
     C        F11.2,' *')
 
316
560   FORMAT (1X,'*',A16,' *   ',I6,'    *   ',I7,'    *',F11.2,' *')
 
317
570   FORMAT (1X,'*',13X,'Wavelength Calibration Diagnosis',14X,'*')
 
318
C Least square solutions for the methods PAIR and ANGLE
 
319
 
 
320
C   Optimal Rotation Angle
 
321
 
 
322
      DET = ABS(N*SYML-SML*SY)
 
323
 
 
324
      IF (DET .GT. 0.01) THEN
 
325
 
 
326
      TGA1 = (N*(SX2*SYML-SXML*SXY)-SX**2*SYML+SML*(SX*SXY-SX2
 
327
     1   *SY)+SX*SXML*SY)/(BIN*(N*(SXML*SY2-SXY*SYML)+SX*SY*SYML+SML*(SX
 
328
     2   Y*SY-SX*SY2)-SXML*SY**2))
 
329
 
 
330
      TGA2 = -(N*SXML-SML*SX)/(BIN*(N*SYML-SML*SY))
 
331
 
 
332
      AOPT = ATAN(TGA1)
 
333
 
 
334
      ELSE
 
335
 
 
336
      AOPT = -999
 
337
   
 
338
      ENDIF
 
339
 
 
340
C   Optimal Binning Factor
 
341
 
 
342
      IF (ABS(A) .GT. 0.02) THEN
 
343
 
 
344
      B1 = COS(A)*(N*(SX2*SYML-SXML*SXY)-SX**2*SYML+SML*(SX*SXY-SX2*SY
 
345
     1   )+SX*SXML*SY)/(SIN(A)*(N*(SXML*SY2-SXY*SYML)+SX*SY*SYML+SML*(SX
 
346
     2   Y*SY-SX*SY2)-SXML*SY**2))
 
347
 
 
348
      ELSE
 
349
 
 
350
      B1 = -999
 
351
 
 
352
      ENDIF
 
353
 
 
354
C   Optimal Absolute Order Number for the first identified line
 
355
 
 
356
      MDP = -((SIN(A)**2*BIN**2*N*SYL-SIN(A)**2*BIN**2*SL*SY+COS(
 
357
     1   A)*SIN(A)*BIN*N*SXL-COS(A)*SIN(A)*BIN*SL*SX)*SYPL+(-SIN(A)**2*B
 
358
     2   IN**2*SPL*SY+COS(A)*SIN(A)*BIN*N*SXPL-COS(A)*SIN(A)*BIN*SPL*SX)
 
359
     3   *SYL+(SIN(A)**2*BIN**2*SL*SPL-SIN(A)**2*BIN**2*N*SPL2)*SY2+SIN(
 
360
     4   A)**2*BIN**2*SPL2*SY**2+(-COS(A)*SIN(A)*BIN*SL*SXPL-COS(A)*SIN(
 
361
     5   A)*BIN*SPL*SXL+2*COS(A)*SIN(A)*BIN*SPL2*SX)*SY+(2*COS(A)*SIN(A)
 
362
     6   *BIN*SL*SPL-2*COS(A)*SIN(A)*BIN*N*SPL2)*SXY+(COS(A)**2*N*SXL-CO
 
363
     7   S(A)**2*SL*SX)*SXPL-COS(A)**2*SPL*SX*SXL+(COS(A)**2*SL*SPL-COS(
 
364
     8   A)**2*N*SPL2)*SX2+COS(A)**2*SPL2*SX**2)/(SIN(A)**2*BIN**2*N*SYL
 
365
     9   **2+(-2*SIN(A)**2*BIN**2*SL*SY+2*COS(A)*SIN(A)*BIN*N*SXL-2*COS(
 
366
     :   A)*SIN(A)*BIN*SL*SX)*SYL+(SIN(A)**2*BIN**2*SL**2-SIN(A)**2*BIN*
 
367
     ;   *2*N*SL2)*SY2+SIN(A)**2*BIN**2*SL2*SY**2+(2*COS(A)*SIN(A)*BIN*S
 
368
     <   L2*SX-2*COS(A)*SIN(A)*BIN*SL*SXL)*SY+(2*COS(A)*SIN(A)*BIN*SL**2
 
369
     =   -2*COS(A)*SIN(A)*BIN*N*SL2)*SXY+COS(A)**2*N*SXL**2-2*COS(A)**2*
 
370
     >   SL*SX*SXL+(COS(A)**2*SL**2-COS(A)**2*N*SL2)*SX2+COS(A)**2*SL2*S
 
371
     ?   X**2)
 
372
 
 
373
      M2OPT = NINT(MDP)
 
374
 
 
375
C Compute new rms if a value is changed.
 
376
 
 
377
      MOLD = M         ! Keep the old value in memory
 
378
      M = M2OPT        ! Replace current value by optimum
 
379
      MODE = 'M'       ! Set the mode
 
380
      GOTO 1000        ! Compute rms 
 
381
1001  CONTINUE         ! Return point
 
382
      M = MOLD         ! Restore old value
 
383
      MRMS = RMSM      ! Store rms
 
384
      MINRMS(1) = RMSM
 
385
 
 
386
      WRITE(TEXT,560) 'Order number',M,M2OPT,MRMS
 
387
      CALL STTPUT(TEXT,ISTAT)
 
388
      CALL STKWRI('OUTPUTI',M2OPT,1,1,KUN,ISTAT)
 
389
 
 
390
      AOLD = A
 
391
      A = AOPT
 
392
      MODE = 'A'
 
393
      GOTO 1000
 
394
1002  CONTINUE         ! Return point
 
395
      A = AOLD
 
396
      ARMS = RMSM
 
397
      MINRMS(2) = RMSM
 
398
 
 
399
      WRITE(TEXT,540) 'Rotation angle',(A*180./PI),(AOPT*180/PI),ARMS
 
400
      CALL STTPUT(TEXT,ISTAT)
 
401
      OUTANG = AOPT*180./PI
 
402
      CALL STKWRR('OUTPUTR',OUTANG,1,1,KUN,ISTAT)
 
403
 
 
404
      BOLD = BIN
 
405
      BIN = B1
 
406
      MODE = 'B'
 
407
      GOTO 1000
 
408
1003  CONTINUE         ! Return point
 
409
      BIN = BOLD
 
410
      BRMS = RMSM
 
411
      MINRMS(3) = RMSM
 
412
 
 
413
      IF (B1 .GT. 0) THEN
 
414
        WRITE(TEXT,540) 'Binning factor',BIN,SQRT(B1),BRMS
 
415
        CALL STTPUT(TEXT,ISTAT)
 
416
      ELSE
 
417
        WRITE(TEXT,540) 'Binning factor',BIN,BIN,0.
 
418
        CALL STTPUT(TEXT,ISTAT)
 
419
      ENDIF
 
420
 
 
421
C
 
422
C Least square solutions for the method PAIR
 
423
 
 
424
      IF (METHOD(1:1) .EQ. 'P') THEN
 
425
 
 
426
      X1  = IX(1)
 
427
      XP1 = IX(2)
 
428
      X2  = IX(3)
 
429
      XP2 = IX(4)
 
430
      P   = IP(3) - IP(1)
 
431
      L1  = IL(1)
 
432
      L2  = IL(3)
 
433
 
 
434
      L1OPT = (((2*L2*M+L2)*P+2*L2*M**2+L2*M)*XP2**2+((L2*P+2*L2*M+2*L2)
 
435
     1   *XP1+((-4*L2*M-2*L2)*P-4*L2*M**2-4*L2*M-L2)*X2+(-L2*P-L2)*X1)
 
436
     2   *XP2+(2*L2*M*P+2*L2*M**2+L2*M)*XP1**2+((L2*P-L2)*X2+
 
437
     3   ((-4*L2*M-2*L2)*P-4*L2*M**2-4*L2*M-L2)*X1)*XP1+((2*L2*M+L2)*P+
 
438
     4   2*L2*M**2+3*L2*M+L2)*X2**2+(-L2*P-2*L2*M)*X1*X2+((2*L2*M+2*L2)
 
439
     5   *P+2*L2*M**2+3*L2*M+L2)*X1**2)/((2*M**2+2*M+2)*XP2**2+(2*M*XP1
 
440
     6   +(-4*M**2-4*M-2)*X2+(-2*M-2)*X1)*XP2+2*M**2*XP1**2+(2*M*X2+(-4
 
441
     7   *M**2-4*M)*X1)*XP1+(2*M**2+2*M+2)*X2**2+(-2*M-2)*X1*X2+(2*M**2
 
442
     8   +4*M+2)*X1**2)
 
443
 
 
444
         L1OLD = L1
 
445
         L1 = L1OPT
 
446
 
 
447
      RMSL1 = -4*((2*L2*P+2*(L2+L1)*M+L2+L1)*(XP2**2+XP1**2+X2**2+X1**2)
 
448
     1 -(XP2+XP1+X2+X1)*(L2*(P+M+1)*XP2+L1*(M+1)*XP1+L2*(P+M)*X2+L1*M*X1
 
449
     2   ))**2/(4*(XP2**2+XP1**2+X2**2+X1**2)-(XP2+XP1+X2+X1)**2)**2-2*(
 
450
     3   XP2+XP1+X2+X1)*(4*(L2*(P+M+1)*XP2+L1*(M+1)*XP1+L2*(P+M)*X2+L1*M
 
451
     4   *X1)-(2*L2*P+2*(L2+L1)*M+L2+L1)*(XP2+XP1+X2+X1))*((2*L2*P+2*(L2
 
452
     5   +L1)*M+L2+L1)*(XP2**2+XP1**2+X2**2+X1**2)-(XP2+XP1+X2+X1)*(L2*(
 
453
     6   P+M+1)*XP2+L1*(M+1)*XP1+L2*(P+M)*X2+L1*M*X1))/(4*(XP2**2+XP1**2
 
454
     7   +X2**2+X1**2)-(XP2+XP1+X2+X1)**2)**2-(4*(L2*(P+M+1)*XP2+L1*(M+1
 
455
     8   )*XP1+L2*(P+M)*X2+L1*M*X1)-(2*L2*P+2*(L2+L1)*M+L2+L1)*(XP2+XP1+
 
456
     9   X2+X1))**2*(XP2**2+XP1**2+X2**2+X1**2)/(4*(XP2**2+XP1**2+X2**2+
 
457
     :   X1**2)-(XP2+XP1+X2+X1)**2)**2+L2**2*(P+M+1)**2+L2**2*(P+M)**2+L
 
458
     ;   1**2*(M+1)**2+L1**2*M**2
 
459
 
 
460
         L1 = L1OLD
 
461
 
 
462
      L2OPT = (((2*L1*M+L1)*P+2*L1*M**2+L1*M)*XP2**2+((L1*P+2*L1*M+2
 
463
     1   *L1)*XP1+((-4*L1*M-2*L1)*P-4*L1*M**2-4*L1*M-L1)*X2+(-L1*P-L1)*X
 
464
     2   1)*XP2+(2*L1*M*P+2*L1*M**2+L1*M)*XP1**2+((L1*P-L1)*X2+((-4*L1*M
 
465
     3   -2*L1)*P-4*L1*M**2-4*L1*M-L1)*X1)*XP1+((2*L1*M+L1)*P+2*L1*M**2+
 
466
     4   3*L1*M+L1)*X2**2+(-L1*P-2*L1*M)*X1*X2+((2*L1*M+2*L1)*P+2*L1*M**
 
467
     5   2+3*L1*M+L1)*X1**2)/((2*P**2+4*M*P+2*M**2)*XP2**2+((2*P+2*M)*XP
 
468
     6   1+(-4*P**2+(-8*M-4)*P-4*M**2-4*M)*X2+(2*P+2*M)*X1)*XP2+(2*P**2+
 
469
     7   (4*M+2)*P+2*M**2+2*M+2)*XP1**2+((-2*P-2*M-2)*X2+(-4*P**2+(-8*M-
 
470
     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**
 
471
     9   2+(-2*P-2*M-2)*X1*X2+(2*P**2+(4*M+2)*P+2*M**2+2*M+2)*X1**2)
 
472
 
 
473
         L2OLD = L2
 
474
         L2 = L2OPT
 
475
 
 
476
      RMSL2 = -4*((2*L2*P+2*(L2+L1)*M+L2+L1)*(XP2**2+XP1**2+X2**2+X1**2)
 
477
     1 -(XP2+XP1+X2+X1)*(L2*(P+M+1)*XP2+L1*(M+1)*XP1+L2*(P+M)*X2+L1*M*X1
 
478
     2   ))**2/(4*(XP2**2+XP1**2+X2**2+X1**2)-(XP2+XP1+X2+X1)**2)**2-2*(
 
479
     3   XP2+XP1+X2+X1)*(4*(L2*(P+M+1)*XP2+L1*(M+1)*XP1+L2*(P+M)*X2+L1*M
 
480
     4   *X1)-(2*L2*P+2*(L2+L1)*M+L2+L1)*(XP2+XP1+X2+X1))*((2*L2*P+2*(L2
 
481
     5   +L1)*M+L2+L1)*(XP2**2+XP1**2+X2**2+X1**2)-(XP2+XP1+X2+X1)*(L2*(
 
482
     6   P+M+1)*XP2+L1*(M+1)*XP1+L2*(P+M)*X2+L1*M*X1))/(4*(XP2**2+XP1**2
 
483
     7   +X2**2+X1**2)-(XP2+XP1+X2+X1)**2)**2-(4*(L2*(P+M+1)*XP2+L1*(M+1
 
484
     8   )*XP1+L2*(P+M)*X2+L1*M*X1)-(2*L2*P+2*(L2+L1)*M+L2+L1)*(XP2+XP1+
 
485
     9   X2+X1))**2*(XP2**2+XP1**2+X2**2+X1**2)/(4*(XP2**2+XP1**2+X2**2+
 
486
     :   X1**2)-(XP2+XP1+X2+X1)**2)**2+L2**2*(P+M+1)**2+L2**2*(P+M)**2+L
 
487
     ;   1**2*(M+1)**2+L1**2*M**2
 
488
 
 
489
 
 
490
         RMSL1 = SQRT(RMSL1/4.)
 
491
         RMSL2 = SQRT(RMSL2/4.)
 
492
 
 
493
        WRITE(TEXT,540) 'Ident. 1 and 2',L1OLD,L1OPT,RMSL1
 
494
        CALL STTPUT(TEXT,ISTAT)
 
495
        WRITE(TEXT,540) 'Ident. 3 and 4',L2OLD,L2OPT,RMSL2
 
496
        CALL STTPUT(TEXT,ISTAT)
 
497
 
 
498
        MINRMS(4) = RMSL1
 
499
        MINRMS(5) = RMSL1
 
500
        MINRMS(6) = RMSL2
 
501
        MINRMS(7) = RMSL2
 
502
        LKOPT(1)  = L1OPT
 
503
        LKOPT(2)  = L1OPT
 
504
        LKOPT(3)  = L2OPT
 
505
        LKOPT(4)  = L2OPT
 
506
 
 
507
      ENDIF
 
508
 
 
509
 
 
510
C  Optimal Wavelength for each of the computed wavelength IL(LOOP)
 
511
 
 
512
      IF (METHOD(1:1) .EQ. 'A') THEN
 
513
 
 
514
 
 
515
      DO 50 K = 1,N
 
516
 
 
517
        PK  = IP(K)
 
518
        LK  = IL(K)
 
519
        XK  = IX(K)
 
520
        YK  = IY(K)
 
521
 
 
522
        SYLK   = 0.D0
 
523
        SYPLK  = 0.D0
 
524
        SPLK   = 0.D0
 
525
        SLK    = 0.D0
 
526
        SL2K   = 0.D0
 
527
        SPL2K  = 0.D0
 
528
        SP2L2K = 0.D0
 
529
        SXLK   = 0.D0
 
530
        SXPLK  = 0.D0
 
531
 
 
532
        DO 30 LOOP = 1 , N
 
533
 
 
534
         IF (LOOP.NE.K) THEN
 
535
           SYLK   = SYLK   + IY(LOOP)*IL(LOOP)
 
536
           SYPLK  = SYPLK  + IY(LOOP)*IP(LOOP)*IL(LOOP)
 
537
           SPLK   = SPLK   + IP(LOOP)*IL(LOOP)
 
538
           SLK    = SLK    + IL(LOOP)
 
539
           SL2K   = SL2K   + IL(LOOP)*IL(LOOP)
 
540
           SPL2K  = SPL2K  + IP(LOOP)*IL(LOOP)**2
 
541
           SP2L2K = SP2L2K + IP(LOOP)**2*IL(LOOP)**2
 
542
           SXLK   = SXLK   + IX(LOOP)*IL(LOOP)
 
543
           SXPLK  = SXPLK  + IX(LOOP)*IP(LOOP)*IL(LOOP)
 
544
         ENDIF
 
545
 
 
546
30      CONTINUE
 
547
 
 
548
C The formula for LKOPT is splitted in two parts (too big for the poor
 
549
C  compiler !!)
 
550
 
 
551
        SBIN = (
 
552
     <   (SIN(A)**2*BIN**2*N*PK+SIN(A)**2*BIN**2*M*N)*YK**2+((2*COS(A)*S
 
553
     =   IN(A)*BIN*N*PK+2*COS(A)*SIN(A)*BIN*M*N)*XK+(-2*SIN(A)**2*BIN**2
 
554
     >   *PK-2*SIN(A)**2*BIN**2*M)*SY+(-2*COS(A)*SIN(A)*BIN*PK-2*COS(A)*
 
555
     ?   SIN(A)*BIN*M)*SX)*YK+(COS(A)**2*N*PK+COS(A)**2*M*N)*XK**2+((-2*
 
556
     @   COS(A)*SIN(A)*BIN*PK-2*COS(A)*SIN(A)*BIN*M)*SY+(-2*COS(A)**2*PK
 
557
     1   -2*COS(A)**2*M)*SX)*XK+((SIN(A)**2*BIN**2-SIN(A)**2*BIN**2*N)*P
 
558
     2   K-SIN(A)**2*BIN**2*M*N+SIN(A)**2*BIN**2*M)*SY2+(SIN(A)**2*BIN**
 
559
     3   2*PK+SIN(A)**2*BIN**2*M)*SY**2+(2*COS(A)*SIN(A)*BIN*PK+2*COS(A)
 
560
     4   *SIN(A)*BIN*M)*SX*SY+((2*COS(A)*SIN(A)*BIN-2*COS(A)*SIN(A)*BIN*
 
561
     5   N)*PK-2*COS(A)*SIN(A)*BIN*M*N+2*COS(A)*SIN(A)*BIN*M)*SXY+((COS(
 
562
     6   A)**2-COS(A)**2*N)*PK-COS(A)**2*M*N+COS(A)**2*M)*SX2+(COS(A)**2
 
563
     7   *PK+COS(A)**2*M)*SX**2)
 
564
 
 
565
       LKOPT(K) = -((SIN(A)**2*BIN**2*N*SYPLK+SIN(A)**2*BIN**2*M*N*SY
 
566
     1   LK+(-SIN(A)**2*BIN**2*SPLK-SIN(A)**2*BIN**2*M*SLK)*SY+COS(A)*SI
 
567
     2   N(A)*BIN*N*SXPLK+COS(A)*SIN(A)*BIN*M*N*SXLK+(-COS(A)*SIN(A)*BIN
 
568
     3   *SPLK-COS(A)*SIN(A)*BIN*M*SLK)*SX)*YK+(COS(A)*SIN(A)*BIN*N*SYPL
 
569
     4   K+COS(A)*SIN(A)*BIN*M*N*SYLK+(-COS(A)*SIN(A)*BIN*SPLK-COS(A)*SI
 
570
     5   N(A)*BIN*M*SLK)*SY+COS(A)**2*N*SXPLK+COS(A)**2*M*N*SXLK+(-COS(A
 
571
     6   )**2*SPLK-COS(A)**2*M*SLK)*SX)*XK+(-SIN(A)**2*BIN**2*SY-COS(A)*
 
572
     7   SIN(A)*BIN*SX)*SYPLK+(-SIN(A)**2*BIN**2*M*SY-COS(A)*SIN(A)*BIN*
 
573
     8   M*SX)*SYLK+(SIN(A)**2*BIN**2*SPLK+SIN(A)**2*BIN**2*M*SLK)*SY2+(
 
574
     9   -COS(A)*SIN(A)*BIN*SXPLK-COS(A)*SIN(A)*BIN*M*SXLK)*SY+(2*COS(A)
 
575
     :   *SIN(A)*BIN*SPLK+2*COS(A)*SIN(A)*BIN*M*SLK)*SXY-COS(A)**2*SX*SX
 
576
     ;   PLK-COS(A)**2*M*SX*SXLK+(COS(A)**2*SPLK+COS(A)**2*M*SLK)*SX2)/
 
577
     +   SBIN
 
578
 
 
579
         LKO = LKOPT(K)
 
580
 
 
581
      RMSLK = -N*((SPLK+M*SLK+LKO*(PK+M))*(SIN(A)**2*BIN**2*SY2+2*COS(A)
 
582
     1   *SIN(A)*BIN*SXY+COS(A)**2*SX2)-(SIN(A)*BIN*SY+COS(A)*SX)*(SIN(A
 
583
     2   )*BIN*(M*(LKO*YK+SYLK)+LKO*PK*YK+SYPLK)+COS(A)*(M*(LKO*XK+SXLK)
 
584
     C   +LKO
 
585
     3   *PK*XK+SXPLK)))**2/(N*(SIN(A)**2*BIN**2*SY2+2*COS(A)*SIN(A)*BIN
 
586
     4   *SXY+COS(A)**2*SX2)-(SIN(A)*BIN*SY+COS(A)*SX)**2)**2-2*(SIN(A)*
 
587
     5 BIN*SY+COS(A)*SX)*(N*(SIN(A)*BIN*(M*(LKO*YK+SYLK)+LKO*PK*YK+SYPLK
 
588
     6)+COS(A)*(M*(LKO*XK+SXLK)+LKO*PK*XK+SXPLK))-(SPLK+M*SLK+LKO*(PK+M)
 
589
     7  )*(SIN(A)*BIN*SY+COS(A)*SX))*((SPLK+M*SLK+LKO*(PK+M))*(SIN(A)**2
 
590
     8   *BIN**2*SY2+2*COS(A)*SIN(A)*BIN*SXY+COS(A)**2*SX2)-(SIN(A)*BIN*
 
591
     9 SY+COS(A)*SX)*(SIN(A)*BIN*(M*(LKO*YK+SYLK)+LKO*PK*YK+SYPLK)+COS(A
 
592
     : )*(M*(LKO*XK+SXLK)+LKO*PK*XK+SXPLK)))/(N*(SIN(A)**2*BIN**2*SY2+2*
 
593
     ;   COS(A)*SIN(A)*BIN*SXY+COS(A)**2*SX2)-(SIN(A)*BIN*SY+COS(A)*SX)*
 
594
     <   *2)**2-(SIN(A)**2*BIN**2*SY2+2*COS(A)*SIN(A)*BIN*SXY+COS(A)**2*
 
595
     + SX2)*(N*(SIN(A)*BIN*(M*(LKO*YK+SYLK)+LKO*PK*YK+SYPLK)+COS(A)*(M*(
 
596
     >LKO*XK+SXLK)+LKO*PK*XK+SXPLK))-(SPLK+M*SLK+LKO*(PK+M))*(SIN(A)*BIN
 
597
     ?   *SY+COS(A)*SX))**2/(N*(SIN(A)**2*BIN**2*SY2+2*COS(A)*SIN(A)*BIN
 
598
     @   *SXY+COS(A)**2*SX2)-(SIN(A)*BIN*SY+COS(A)*SX)**2)**2+2*M*SPL2K+
 
599
     1   SP2L2K+M**2*SL2K+LKO**2*(PK+M)**2
 
600
 
 
601
 
 
602
         RMSLK = SQRT(RMSLK/N)
 
603
         MINRMS(3+K) = RMSLK
 
604
         WRITE(TEXT,550) K,LK,LKOPT(K),RMSLK
 
605
         CALL STTPUT(TEXT,ISTAT)
 
606
 
 
607
50      CONTINUE
 
608
 
 
609
        ENDIF
 
610
 
 
611
        WRITE(TEXT,510)
 
612
        CALL STTPUT(TEXT,ISTAT)
 
613
 
 
614
        GOTO 2000
 
615
 
 
616
C **********************************************************************
 
617
 
 
618
1000    CONTINUE
 
619
C       Computes RMS(A,M,BIN). There are too many parameters to create
 
620
C       a subroutine. Let's use the good old GOTO.
 
621
 
 
622
       RMSM = -N*((SPL+M*SL)*(SIN(A)**2*BIN**2*SY2+2*COS(A)*SIN(A)*BIN*S
 
623
     1   XY+COS(A)**2*SX2)-(SIN(A)*BIN*SY+COS(A)*SX)*(SIN(A)*BIN*(SYPL+M
 
624
     2   *SYL)+COS(A)*(SXPL+M*SXL)))**2/(N*(SIN(A)**2*BIN**2*SY2+2*COS(A
 
625
     3   )*SIN(A)*BIN*SXY+COS(A)**2*SX2)-(SIN(A)*BIN*SY+COS(A)*SX)**2)**
 
626
     4   2-2*(SIN(A)*BIN*SY+COS(A)*SX)*(N*(SIN(A)*BIN*(SYPL+M*SYL)+COS(A
 
627
     5   )*(SXPL+M*SXL))-(SPL+M*SL)*(SIN(A)*BIN*SY+COS(A)*SX))*((SPL+M*S
 
628
     6   L)*(SIN(A)**2*BIN**2*SY2+2*COS(A)*SIN(A)*BIN*SXY+COS(A)**2*SX2)
 
629
     7   -(SIN(A)*BIN*SY+COS(A)*SX)*(SIN(A)*BIN*(SYPL+M*SYL)+COS(A)*(SXP
 
630
     8   L+M*SXL)))/(N*(SIN(A)**2*BIN**2*SY2+2*COS(A)*SIN(A)*BIN*SXY+COS
 
631
     9   (A)**2*SX2)-(SIN(A)*BIN*SY+COS(A)*SX)**2)**2-(SIN(A)**2*BIN**2*
 
632
     :   SY2+2*COS(A)*SIN(A)*BIN*SXY+COS(A)**2*SX2)*(N*(SIN(A)*BIN*(SYPL
 
633
     ;   +M*SYL)+COS(A)*(SXPL+M*SXL))-(SPL+M*SL)*(SIN(A)*BIN*SY+COS(A)*S
 
634
     <   X))**2/(N*(SIN(A)**2*BIN**2*SY2+2*COS(A)*SIN(A)*BIN*SXY+COS(A)*
 
635
     =   *2*SX2)-(SIN(A)*BIN*SY+COS(A)*SX)**2)**2+2*M*SPL2+SP2L2+M**2*SL
 
636
     >   2
 
637
 
 
638
         RMSM = SQRT(RMSM/N)
 
639
 
 
640
        IF (MODE(1:1) .EQ. 'M') GOTO 1001  ! Absolute order number
 
641
        IF (MODE(1:1) .EQ. 'A') GOTO 1002  ! Angle
 
642
        IF (MODE(1:1) .EQ. 'B') GOTO 1003  ! Bin
 
643
 
 
644
C **********************************************************************
 
645
 
 
646
1900    CONTINUE
 
647
        ERROR = 1
 
648
 
 
649
2000    CONTINUE
 
650
        CALL STKWRI('OUTPUTI',ERROR,1,1,KUN,ISTAT)
 
651
        CALL STSEPI
 
652
        STOP
 
653
        END
 
654
 
 
655
 
 
656
 
 
657
 
 
658