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

« back to all changes in this revision

Viewing changes to stdred/echelle/proc/neciden.prg

  • 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:
 
1
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++
 
2
!.COPYRIGHT   (C) 1991-2011 European Southern Observatory
 
3
!.IDENT       neciden.prg
 
4
!.AUTHOR      Pascal Ballester,  ESO - Garching
 
5
!.KEYWORDS    Spectroscopy, Echelle, 
 
6
!.PURPOSE     Wavelength calibration of echelle spectra
 
7
!             Implement command IDENT/ECHELLE
 
8
!    IDENT/ECH  WLC  Lincat  Degree  Tol  Initol,Dev,Stop  Starter Session
 
9
!.VERSION     1.0    Creation    18.01.91
 
10
! 110525        last modif
 
11
 
12
!-------------------------------------------------------
 
13
!
 
14
IF P1(1:1) .NE. "?"   WRITE/KEY    WLC/C/1/60      {P1}
 
15
IF P2(1:1) .NE. "?"   WRITE/KEY    LINCAT/C/1/60   {P2}
 
16
IF P3(1:1) .NE. "?"   DC = {P3}
 
17
IF P4(1:1) .NE. "?"   TOL = {P4}
 
18
IF P5(1:1) .NE. "?"   WRITE/KEY    WLCITER/R/1/5   {P5}
 
19
IF P6(1:1) .NE. "?"   WRITE/KEY    WLCMTD/C/1/10   {P6}
 
20
IF P7(1:1) .NE. "?"   WRITE/KEY    GUESS/C/1/60    {P7}
 
21
IF P8(1:1) .NE. "?"   WRITE/KEY    CCDBIN/R/1/2    {P8}
 
22
 
 
23
DEFINE/LOCAL    METHOD/C/1/7  "YYYYYYY" ?  +lower
 
24
IF WLCMTD(1:1) .EQ.  "P"   WRITE/KEY  METHOD/C/1/7  "PAIR   "
 
25
IF WLCMTD(1:1) .EQ.  "A"   WRITE/KEY  METHOD/C/1/7  "ANGLE  "
 
26
IF WLCMTD(1:1) .EQ.  "T"   WRITE/KEY  METHOD/C/1/7  "TWO-D  "
 
27
IF WLCMTD(1:1) .EQ.  "G"   WRITE/KEY  METHOD/C/1/7  "GUESS  "
 
28
IF WLCMTD(1:1) .EQ.  "R"   WRITE/KEY  METHOD/C/1/7  "RESTART"
 
29
IF WLCMTD(1:1) .EQ.  "O"   WRITE/KEY  METHOD/C/1/7  "ORDER  "
 
30
IF WLCMTD(1:1) .EQ.  "L"   WRITE/KEY  METHOD/C/1/7  "LINEAR "
 
31
 
 
32
IF METHOD(1:1) .EQ. "Y" THEN
 
33
   WRITE/OUT "IDENT/ECHELLE: unknown method {WLCMTD}"
 
34
   WRITE/OUT "Method can be PAIR, ANGLE, TWO-D, GUESS, RESTART, ORDER, LINEAR"
 
35
   RETURN/EXIT
 
36
ENDIF
 
37
 
 
38
VERIFY/ECHELLE {WLC}
 
39
VERIFY/ECHELLE {LINTAB} TAB
 
40
VERIFY/SPEC    {LINCAT} MID_ARC LINCAT TAB
 
41
 
 
42
IF WLCVISU(1:1) .EQ. "Y"   THEN
 
43
  DISPLAY/ECHELLE {WLC}
 
44
  CLEAR/CHAN OVER
 
45
ENDIF
 
46
 
 
47
WRITE/KEY  IN_A {LINTAB}
 
48
WRITE/KEY  IN_B {LINCAT}
 
49
INPUTR(1) = -10 ! Verification mode
 
50
RUN STD_EXE:spematch 
 
51
IF INPUTR(1) .EQ. -1    THEN
 
52
    WRITE/OUT   "Now sorting table {LINCAT}."
 
53
    SORT/TABLE  {LINCAT}  :WAVE
 
54
    WRITE/KEY  IN_A {LINTAB}
 
55
    WRITE/KEY  IN_B {LINCAT}
 
56
    INPUTR(1) = -10 ! Verification mode
 
57
    RUN STD_EXE:spematch 
 
58
ENDIF
 
59
IF INPUTR(1) .LT. 0  RETURN/EXIT
 
60
 
 
61
! Initial check : Number of elements of LINE > 2 * NBORD * (DEGREE + 1)
 
62
IF DC .GT.  8   ERROR/ECHELLE  IDENT/ECHELLE   DC
 
63
IF WLCOPT(1:1) .NE. "1" .AND. WLCOPT(1:1) .NE. "2" THEN
 
64
    ERROR/ECHELLE  IDENT/ECHELLE   WLCOPT
 
65
ENDIF
 
66
IF WLCREG(1:1) .NE. "S" .AND. WLCREG(1:1) .NE. "R" then
 
67
   write/out "Error: Keyword WLCREG = {WLCREG}. Accepted values: Standard, Robust"
 
68
   error/echelle  IDENT/ECHELLE WLCREG
 
69
ENDIF
 
70
 
 
71
DEFINE/LOCAL   WLCMTD/C/1/12 {WLCMTD} ? +lower 
 
72
DEFINE/LOCAL   DEGREE/I/1/1  {DC}     ? +lower  
 
73
! Degree of poly. for orders fit.
 
74
DEFINE/LOCAL   NBORD/I/1/1   {ECHORD(1)} ? +lower   ! Number of orders
 
75
 
 
76
DEFINE/LOCAL ROTNAM/C/1/6  ROTATE     ? +lower 
 
77
DEFINE/LOCAL   nbr/r/1/3     0,0,0    ? +lower 
 
78
DEFINE/LOCAL   last/i/1/1    0        ? +lower 
 
79
DEFINE/LOCAL   line/i/1/1    0        ? +lower 
 
80
DEFINE/LOCAL   xymax/r/1/2   0,0      ? +lower 
 
81
DEFINE/LOCAL   nbsel/i/1/1   0        ? +lower 
 
82
DEFINE/LOCAL   alpha/d/1/3   0,0,0    ? +lower 
 
83
DEFINE/LOCAL   coupl/i/1/1   0        ? +lower 
 
84
DEFINE/LOCAL   ordiff/i/1/1  0        ? +lower  
 
85
DEFINE/LOCAL   ordpos/i/1/2  0,0      ? +lower 
 
86
DEFINE/LOCAL   rms/r/1/2     0,0      ? +lower 
 
87
DEFINE/LOCAL   pixel/r/1/3   0,0,0    ? +lower 
 
88
DEFINE/LOCAL   select/r/1/3  0,0,0    ? +lower 
 
89
DEFINE/LOCAL   alpmax/r/1/1  0        ? +lower 
 
90
DEFINE/LOCAL   RATIO/R/1/1   0.       ? +lower 
 
91
DEFINE/LOCAL   NBRMIN/I/1/1  0        ? +lower 
 
92
DEFINE/LOCAL   NBLMIN/I/1/1  0        ? +lower 
 
93
DEFINE/LOCAL   LWCL/R/1/5    0.,0.,0.,0.,0. ? +lower
 
94
define/local   wlcstep/D/1/1 0.       ? +lower
 
95
 
 
96
WRITE/KEYWORD  xypos/r/1/2   0,0
 
97
WRITE/KEYWORD  wave/d/1/1    0
 
98
WRITE/KEYWORD  sequ/I/1/1    0
 
99
WRITE/KEYWORD  ordval/r/1/1  0
 
100
WRITE/KEYWORD  answ/c/1/3    HAI
 
101
COPY/KEYWORD   wlciter/R/1/5 LWCL/R/1/5
 
102
copy/dk        {wlc} step/D/1/1  wlcstep
 
103
 
 
104
SET/FORMAT  I1  F12.5
 
105
 
 
106
WRITE/OUT      "*** Wavelength calibration ***"            
 
107
WRITE/OUT      "   "                                        
 
108
 
 
109
IF OUTMODE(1:1) .NE. "S" THEN ! Mode Silent
 
110
WRITE/OUT      " Method                    : {WLCMTD}"      
 
111
WRITE/OUT      " Dispersion relation type  : {WLCOPT}"
 
112
WRITE/OUT      " Guess session (if any)    : {GUESS}"       
 
113
WRITE/OUT      " WLC frame                 : {WLC}"         
 
114
WRITE/OUT      " Line catalog              : {LINCAT}"      
 
115
WRITE/OUT      " Number of orders          : {NBORD}"       
 
116
WRITE/OUT      " Polynomial degree         : {DEGREE}"      
 
117
WRITE/OUT      " Binning                   : {CCDBIN(1)},{CCDBIN(2)}"    
 
118
WRITE/OUT      "** Identification loop **"    
 
119
WRITE/OUT      " Initial accuracy          : {LWCL(1)} pixel"    
 
120
WRITE/OUT      " Next neighbour parameter  : {LWCL(2)}"          
 
121
WRITE/OUT      " Maximal error             : {LWCL(3)} pixel"
 
122
WRITE/OUT      " Initial error             : {LWCL(4)} pixel"
 
123
WRITE/OUT      " Iteration error           : {LWCL(5)} std dev."
 
124
WRITE/OUT      "** Rejection loop **"                               
 
125
WRITE/OUT      " Tolerance                 : {TOL}"                 
 
126
ENDIF
 
127
 
 
128
IF LWCL(2) .LE. 0.  ERROR/ECHELLE  IDENTIFY/ECHELLE  LWCL(2)
 
129
 
 
130
WRITE/KEY      LAST/I/1/1   {{LINTAB},TBLCONTR(4)} ! Number of elements in LINE
 
131
LINE   = (2. * NBORD * (DEGREE + 2.))  ! Minimum number of elements in table line
 
132
NBRMIN = 2*(DEGREE+1)*(DEGREE+1)       ! Minimum number of identifications to avoid echelle relat.
 
133
 
 
134
IF LAST .LT. LINE  THEN
 
135
         WRITE/OUT "WARNING : The table LINE, produced by SEARCH/ECHELLE"
 
136
         WRITE/OUT "perhaps does not contain enough lines for this calibration."
 
137
         WRITE/OUT "Corrective actions can be to repeat:"
 
138
         WRITE/OUT "       - SEARCH/ECH with a lower threshold"
 
139
         WRITE/OUT "       - IDENT/ECH with a lower degree"
 
140
ENDIF
 
141
 
 
142
IF WLCMTD(1:1) .EQ. "O"  THEN
 
143
   COPY/DK  {LINTAB} PIXEL  PIXEL 
 
144
   GOTO  END
 
145
ENDIF
 
146
 
 
147
WRITE/DESCR     {LINTAB}      HISTORY_UPDA/I/1/1   0    ! No update of history
 
148
 
 
149
CREATE/COLUMN   {LINTAB}      :ORDER       "Abs. order"            I6   {SESSOUTV}
 
150
CREATE/COLUMN   {LINTAB}      :IDENT       "Wavelength"    F10.4   R*8  {SESSOUTV}
 
151
CREATE/COLUMN   {LINTAB}      :WAVEC       "Wavelength"    F10.4   R*8  {SESSOUTV}
 
152
CREATE/COLUMN   {LINTAB}      :AUX         "ORD.lambda"    F10.4   R*8  {SESSOUTV}
 
153
CREATE/COLUMN   {LINTAB}      :XROT        "Rotat. pos"    F7.2         {SESSOUTV}
 
154
 
 
155
CREATE/COLUMN   {LINTAB}      :ERROR       "Wav. un.  "    F10.6   R*8  {SESSOUTV}
 
156
CREATE/COLUMN   {LINTAB}      :DISTANCE    "Wav. un.  "    F10.6   R*8  {SESSOUTV}
 
157
CREATE/COLUMN   {LINTAB}      :RESIDUAL    "Wavel. un."    F10.6   R*8  {SESSOUTV}
 
158
CREATE/COLUMN   {LINTAB}      :SELECT      "Logical   "    I6      I    {SESSOUTV}
 
159
 
 
160
CREATE/TABLE    &t    1    20
 
161
CREATE/COLUMN   &t    :COL
 
162
 
 
163
COMPUTE/TABLE   {LINTAB}      :ERROR = -1.
 
164
COMPUTE/TABLE   {LINTAB}      :IDENT = 0.
 
165
COMPUTE/TABLE   {LINTAB}      :RESIDUAL = 0.
 
166
COMPUTE/TABLE   {LINTAB}      :SEQU  = sequence
 
167
COMPUTE/TABLE   {LINTAB}      :SELECT = 0.
 
168
 
 
169
IF  WLCMTD(1:1) .EQ. "R"  THEN
 
170
    COPY/DK  middummr.tbl  WLCMTD   WLCMTD
 
171
ELSE
 
172
    CREATE/TABLE    &r    4    10
 
173
    CREATE/COLUMN   &r    :X
 
174
    CREATE/COLUMN   &r    :Y
 
175
    CREATE/COLUMN   &r    :DX
 
176
    CREATE/COLUMN   &r    :DY
 
177
ENDIF
 
178
!
 
179
! Use different starters depending on the configuration
 
180
! Interactive PAIR method is supported by the routine PAIR
 
181
! Other interactive methods (ANGLE, TWO-D) are supported by ANGLE
 
182
 
 
183
IF WLCMTD(1:1) .EQ. "P" @s neciden,pair
 
184
 
 
185
IF WLCMTD(1:1) .EQ. "A" @s neciden,angle
 
186
IF WLCMTD(1:1) .EQ. "T" @s neciden,angle
 
187
 
 
188
IF WLCMTD(1:1) .EQ. "G" @s neciden,guess
 
189
IF WLCMTD(1:1) .EQ. "L" @s neciden,hough
 
190
 
 
191
COPY/KD  PIXEL  {LINTAB}     PIXEL
 
192
COPY/KD  WLCMTD middummr.tbl WLCMTD
 
193
 
 
194
! Branch to the main iteration loop
 
195
ITERATION:
 
196
IF WLCREG(1:1) .EQ. "S"  @s neciden,iter
 
197
IF WLCREG(1:1) .EQ. "R"  @s neciden,htiter
 
198
 
 
199
END:
 
200
pixel(1) = M$ABS(pixel(1))
 
201
!
 
202
! Now, to be able to find estimated polynomials where there is not 
 
203
! enough lines, one fits a 2D polynomial as m.lambda=f(x,y). The reason 
 
204
! why this is not used from the beginning is that one needs 
 
205
! (DEGREE+2)*(DEGREE+2) observation values (at least).
 
206
!
 
207
COMPUTE/TABLE   {LINTAB}   :AUX = 0.
 
208
COMPUTE/TABLE   {LINTAB}   :AUX = :ORDER * :IDENT
 
209
SELECT/TABLE    {LINTAB}   (:IDENT.GT.0.0)  {SESSOUTV}
 
210
REGRES/POLY     {LINTAB}   :AUX   :X,:ORDER  {DEGREE},{DEGREE}  KEYLONG {SESSOUTV}
 
211
SAVE/REGR       {LINTAB}   REGR   KEYLONG 
 
212
 
 
213
IF TOL .EQ. 0. THEN
 
214
   PIXEL(2) = 3.*KEYLONGR(5)*2./(echord(2)+echord(3))
 
215
ELSE
 
216
 IF TOL .GT. 0. THEN
 
217
     PIXEL(2)        =          TOL*PIXEL(1)
 
218
 ELSE
 
219
     PIXEL(2)        =          M$ABS(TOL)
 
220
 ENDIF
 
221
ENDIF
 
222
 
 
223
COMPUTE/REGRE {LINTAB} :WAVEC = REGR
 
224
COMPUTE/TABLE {LINTAB} :WAVEC = :WAVEC/:ORDER
 
225
COMPUTE/TABLE {LINTAB} :RESIDUAL = :WAVEC - :IDENT
 
226
SELECT/TABLE  {LINTAB} ABS(:RESIDUAL) .LT. {PIXEL(2)} {SESSOUTV}
 
227
 
 
228
Write/Out     "Final Selection: Threshold = {pixel(2)} wav. units, Nb = {outputi(1)}"
 
229
IF WLCVISU(1:1) .EQ. "Y" THEN
 
230
   CLEAR/CHAN OVER
 
231
   LOAD/TABLE {LINTAB} :X :YNEW
 
232
ENDIF
 
233
 
 
234
! Save the final dispersion relation and order numbers for later
 
235
! use in GUESS mode.
 
236
REGRES/POLY  {LINTAB} :AUX :X,:ORDER  {DEGREE},{DEGREE}  KEYLONG {SESSOUTV}
 
237
SAVE/REGR    {LINTAB} REGR KEYLONG 
 
238
REGRES/POLY  {LINTAB} :ORDER :X,:YNEW  {DEGREE},{DEGREE}  KEYLONG {SESSOUTV}
 
239
SAVE/REGR    {LINTAB} RORD KEYLONG 
 
240
 
 
241
SELECT/TABLE    {LINTAB}   (:IDENT.GT.0.0) {SESSOUTV}
 
242
IF WLCOPT(1:1) .EQ. "2" SELECT/TABLE {LINTAB} SELECT.AND.ABS(:RESIDUAL).LT.{PIXEL(2)}
 
243
 
 
244
INPUTI(1)   = DC
 
245
INPUTI(2)   = echord(1)
 
246
INPUTI(3)   = imsize(1)
 
247
INPUTI(4)   = {WLCOPT(1:1)}
 
248
 
 
249
INPUTR(1)   = PIXEL(2)  ! Tolerance in wavelength units
 
250
 
 
251
INPUTD(1)   = {{WLC},START(1)}
 
252
INPUTD(2)   = {{WLC},STEP(1)}
 
253
 
 
254
RUN  STD_EXE:neciden
 
255
 
 
256
WLRANGE(1) = OUTPUTR(1)
 
257
WLRANGE(2) = OUTPUTR(2)
 
258
 
 
259
@s neciden,verif
 
260
 
 
261
! Estimate value of K for optional IUE-like RIPPLE correction
 
262
set/format F20.10
 
263
STAT/TAB  {LINTAB}   :AUX  {SESSOUTV}
 
264
write/key  SAMPLE  {{LINTAB},PIXEL(1)}
 
265
RIPK = OUTPUTR(3)
 
266
echord(6) = echord(3) + 1
 
267
write/out "Set parameters SAMPLE={SAMPLE} and RIPK={RIPK}"
 
268
 
 
269
EXIT:
 
270
 
 
271
if wlcvisu(1:1) .eq. "Y" PLOT/RESIDUAL
 
272
 
 
273
RETURN
 
274
 
 
275
!*****************************************************************************
 
276
 
 
277
ENTRY PAIR
 
278
 
 
279
INTERACT:
 
280
 
 
281
IF METHOD(1:1) .NE. "R" THEN
 
282
 
 
283
WRITE/OUT    "Info: You must identify lines which are repeated in"
 
284
WRITE/OUT    "      overlapped regions of two adjacent orders. You"
 
285
WRITE/OUT    "      will do it two times."
 
286
WRITE/OUT    "    "
 
287
WRITE/OUT    "With the cursor, select two lines that appear in two "
 
288
WRITE/OUT    "different orders each. Click them in the following order :"
 
289
WRITE/OUT    "   - First line: left-hand side, then right-side"
 
290
WRITE/OUT    "   - Second line: left-hand-side, then right-side"
 
291
 
 
292
ENDIF
 
293
 
 
294
@s neciden,cursor {METHOD}
 
295
@s neciden,overlay
 
296
 
 
297
if outputi(1) .ne. 3 then
 
298
   write/out "you haven't selected four lines. please restart."
 
299
   write/descr  middummr.tbl  tblcontr/i/4/1   0
 
300
   GOTO  INTERACT
 
301
endif
 
302
 
 
303
! *** Find selections in table LINE
 
304
 
 
305
 
 
306
@s neciden,preploop
 
307
 
 
308
DO line = 1  'last'
 
309
 
 
310
     @s neciden,assoc
 
311
 
 
312
     IF line .ne. 2 THEN
 
313
        IF line .ne. 4 THEN
 
314
           COUPL(1)   =    {{LINTAB},:ORDER,@{SEQU}}
 
315
           INQUIRE/KEY     WAVE/D/1/1  "Sequence no. {LINE}, Order no. {COUPL}. Enter wavelength : "
 
316
        ENDIF
 
317
     ENDIF
 
318
     COPY/KT               WAVE/D/1/1  {LINTAB}   :IDENT   @{SEQU}
 
319
 
 
320
ENDDO
 
321
 
 
322
! computes rotation angle and peforms rotation
 
323
ANGLE:
 
324
 
 
325
alpha(1) = ({middummr,:y,@1}-{middummr,:y,@2})/({middummr,:x,@1}-{middummr,:x,@2})
 
326
alpha(2) = ({middummr,:y,@3}-{middummr,:y,@4})/({middummr,:x,@3}-{middummr,:x,@4})
 
327
 
 
328
ALPHA(1) = M$ATAN (ALPHA(1))
 
329
ALPHA(2) = M$ATAN (ALPHA(2))
 
330
 
 
331
inputr(1) = 0.
 
332
! @s echidenew,analys
 
333
@s neciden,analys
 
334
alpmax = outputr(1)
 
335
 
 
336
! Average angle as a result of geometrical analysis
 
337
ALPHA(3) =  (0.5)*(ALPHA(1) + ALPHA(2))
 
338
WRITE/OUT  "      >>>   Geometrical orientation analysis   >>>  Angle: 'ALPHA(3)' deg."
 
339
WRITE/OUT  "      >>>   Analytical  orientation analysis   >>>  Angle: 'ALPMAX' deg."
 
340
 
 
341
define/local angdif/R/1/1 0.
 
342
angdif = M$ABS(ALPHA(3)-ALPMAX)
 
343
IF angdif .gt. 8. then
 
344
    Write/Out "Warning: An inconsistency has been detected between geometrical"
 
345
    Write/Out "         and analytical determination of the rotation angle."
 
346
    Write/Out "         This could result from misidentifications or wrong "
 
347
    Write/Out "         order number. "
 
348
    INQUIRE/KEYWORD ANSW/C/1/3   "Continue anyway (y/n,def=n) ? "
 
349
    IF ANSW(1:1) .EQ. "N" .OR. AUX_MODE(7) .EQ. 0  RETURN/EXIT
 
350
ENDIF
 
351
 
 
352
!Prepares table
 
353
COMPUTE/TABLE  {LINTAB}  :XROT = :X*{CCDBIN(1)}*COS('alpha(3)')+:YNEW*{CCDBIN(2)}*SIN('alpha(3)')
 
354
 
 
355
 
 
356
COMPUTE/TABLE  {LINTAB}  :AUX = :ORDER * :IDENT
 
357
SELECT/TABLE   {LINTAB}  (:IDENT.GT.0.0) {SESSOUTV}
 
358
NBR(1) = OUTPUTI(1)
 
359
REGRES/POLY    {LINTAB}  :AUX   :XROT  2 {SESSOUTV}
 
360
SAVE/REGRES    {LINTAB}  REGR
 
361
 
 
362
RMS(1)   =  2.  *  OUTPUTR(5) / (ORDPOS(1) + ORDPOS(2))
 
363
PIXEL(1) =  2.  *  OUTPUTD(2) / (ORDPOS(1) + ORDPOS(2))
 
364
PIXEL(1) = M$ABS(PIXEL(1)*wlcstep)
 
365
 
 
366
IF {LWCL(1)} .GT. 0. THEN
 
367
        PIXEL(2) =  LWCL(1)*PIXEL(1)
 
368
ELSE
 
369
        PIXEL(2) = M$ABS(LWCL(1))
 
370
ENDIF
 
371
 
 
372
WRITE/OUT  "*** Initial, coarse estimate ***"
 
373
WRITE/OUT  "Mean pixel width      : {PIXEL(1)} wavel. units"
 
374
WRITE/OUT  "RMS error of relation : {RMS(1)} wavel. units"
 
375
WRITE/OUT  "Maximum allowed error : {PIXEL(2)} wavel. units"
 
376
 
 
377
IF rms(1) .gt. pixel(2) THEN
 
378
   WRITE/OUT "This accuracy appears insufficient to attempt automatic line identification."
 
379
   WRITE/OUT "Please check your identifications."
 
380
 
 
381
   INPUTR(1) = ALPHA(3)
 
382
   INPUTR(2) = 1.
 
383
   @s neciden,analys
 
384
   RATIO = M$ABS(RMS(1)/PIXEL(2))
 
385
   IF RATIO .LT. LWCL(3)  THEN
 
386
   WRITE/OUT "However the accuracy may be still acceptable."
 
387
   INQUIRE/KEY   ANSW/C/1/3  "Start anyway (y/n, default no) ?"
 
388
   IF ANSW(1:1) .EQ. "Y" THEN
 
389
           GOTO END
 
390
   ELSE
 
391
           RETURN/EXIT
 
392
           WRITE/DESCR  middummc.tbl  TBLCONTR/I/4/1   0
 
393
           IF METHOD(1:1) .NE. "R"  WRITE/DESCR  middummr.tbl TBLCONTR/I/4/1 0
 
394
           GOTO  INTERACT
 
395
   ENDIF
 
396
   ELSE
 
397
   WRITE/OUT "Inaccuracy is too large. Aborted..."
 
398
   RETURN/EXIT
 
399
   ENDIF
 
400
ENDIF
 
401
 
 
402
END:
 
403
 !Prepares new table with analytical rotation angle
 
404
 COMPUTE/TABLE  {LINTAB}  :XROT = :X*{CCDBIN(1)}*COS('alpmax')+:YNEW*{CCDBIN(2)}*SIN('alpmax')
 
405
 
 
406
 COMPUTE/TABLE  {LINTAB}  :AUX = :ORDER * :IDENT
 
407
 SELECT/TABLE   {LINTAB}  (:IDENT.GT.0.0)
 
408
 NBR(1) = OUTPUTI(1)
 
409
 REGRES/POLY    {LINTAB}  :AUX   :XROT  2
 
410
 SAVE/REGRES    {LINTAB}  REGR
 
411
 
 
412
 RMS(1)   =  2.  *  OUTPUTR(5) / (ORDPOS(1) + ORDPOS(2))
 
413
 PIXEL(1) =  2.  *  OUTPUTD(2) / (ORDPOS(1) + ORDPOS(2))
 
414
 pixel(1) =  M$ABS(pixel(1))
 
415
 
 
416
 IF {LWCL(1)} .GT. 0. THEN
 
417
        PIXEL(2) =  {LWCL(1)} *  PIXEL(1)
 
418
 ELSE
 
419
        PIXEL(2) = M$ABS({LWCL(1)})
 
420
 ENDIF
 
421
 
 
422
RETURN
 
423
 
 
424
! ***************** Start with 4 lines and play with rotation ************
 
425
 
 
426
ENTRY ANGLE
 
427
 
 
428
DEFINE/LOCAL  DEGSTART/I/1/1  2  ? +lower 
 
429
 
 
430
METH2:
 
431
 
 
432
IF WLCMTD(1:1) .EQ. "A"  NBLMIN = 4
 
433
IF WLCMTD(1:1) .EQ. "C"  NBLMIN = DEGSTART + 2
 
434
IF WLCMTD(1:1) .EQ. "T"  NBLMIN = (DEGSTART+1)*(DEGSTART+1)+1
 
435
 
 
436
SET/FORMAT  I1
 
437
 
 
438
IF METHOD(1:1) .NE. "R" THEN
 
439
WRITE/OUT    "Point at least {NBLMIN} lines of which you know the wavelength"
 
440
ENDIF
 
441
 
 
442
@s neciden,cursor {METHOD}
 
443
@s neciden,overlay
 
444
 
 
445
outputi(1) = outputi(1) + 1
 
446
 
 
447
IF outputi(1) .lt. NBLMIN THEN
 
448
   WRITE/OUT "You have selected less than {NBLMIN} lines. Please start again."
 
449
   WRITE/DESCR middummr.tbl  TBLCONTR/I/4/1  0
 
450
   GOTO METH2
 
451
ENDIF
 
452
SET/FORMAT
 
453
 
 
454
@s neciden,preploop
 
455
 
 
456
DO line = 1  'last'
 
457
 
 
458
     @s neciden,assoc
 
459
 
 
460
     COUPL(1)   =    {{LINTAB},:ORDER,@{SEQU}}
 
461
     INQUIRE/KEY     WAVE/D/1/1  "Sequence no. {LINE}, Order no. 'COUPL'. Enter wavelength : "
 
462
     COPY/KT         WAVE/D/1/1  {LINTAB}   :IDENT   @{SEQU}
 
463
 
 
464
ENDDO
 
465
 
 
466
!KB SELECT/TABLE  {LINTAB}  (:IDENT.NE.0.0) {SESSOUTV}
 
467
SELECT/TAB       {LINTAB} (ABS(:IDENT).GT.0.00001)  {SESSOUTV}
 
468
 
 
469
IF WLCMTD(1:1) .EQ. "A" THEN ! Method ANGLE
 
470
  inputr(1) = 0.
 
471
  @s neciden,analys
 
472
  alpha(3) = outputr(1)
 
473
  COMPUTE/TABLE  {LINTAB}  :XROT=:X*COS('ALPHA(3)')+:YNEW*SIN('ALPHA(3)')
 
474
  COMPUTE/TABLE  {LINTAB}  :AUX = :ORDER * :IDENT
 
475
  SELECT/TABLE   {LINTAB}  (:IDENT.GT.0.0)
 
476
  NBR(1) = OUTPUTI(1)
 
477
  REGRES/POLY    {LINTAB}  :AUX   :XROT  {DEGSTART}
 
478
  PIXEL(1) = M$ABS(OUTPUTD(2)*wlcstep)
 
479
ENDIF
 
480
 
 
481
IF WLCMTD(1:1) .EQ. "T" THEN   ! Method TWO-D
 
482
  COMPUTE/TABLE  {LINTAB}  :AUX = :ORDER * :IDENT
 
483
  SELECT/TABLE   {LINTAB}  (:IDENT.GT.0.0)
 
484
  NBR(1) = OUTPUTI(1)
 
485
  REGRES/POLY    {LINTAB}  :AUX  :X,:ORDER  {DEGSTART},{DEGSTART}
 
486
  PIXEL(1) = M$ABS(OUTPUTD(2)*wlcstep)
 
487
ENDIF
 
488
 
 
489
SAVE/REGRES    {LINTAB}  REGR
 
490
RMS(1)   = OUTPUTR(5)
 
491
IF WLCMTD(1:1) .NE. "C"  THEN
 
492
   RMS(1)   =  2.*RMS(1)/(ORDPOS(1) + ORDPOS(2))
 
493
   PIXEL(1) =  2.*PIXEL(1)/ (ORDPOS(1) + ORDPOS(2))
 
494
ENDIF
 
495
 
 
496
IF {LWCL(1)} .GT. 0. THEN
 
497
        PIXEL(2) =  LWCL(1)*  PIXEL(1)
 
498
ELSE
 
499
        PIXEL(2) = M$ABS(LWCL(1))
 
500
ENDIF
 
501
 
 
502
WRITE/OUT  "*** Initial, coarse estimate ***"
 
503
WRITE/OUT  "Mean pixel width      : {PIXEL(1)} wavel. units"
 
504
WRITE/OUT  "RMS error of relation : {RMS(1)} wavel. units"
 
505
WRITE/OUT  "Maximum allowed error : {PIXEL(2)} wavel. units"
 
506
IF rms(1) .gt. pixel(2) THEN
 
507
   WRITE/OUT "This accuracy appears insufficient to attempt automatic line identification."
 
508
   WRITE/OUT "Please check your identifications."
 
509
   IF WLCMTD(1:1) .EQ. "A" THEN
 
510
      inputr(1) = alpha(3)
 
511
      @s neciden,analys
 
512
      alpha(3) = outputr(1)
 
513
   ENDIF
 
514
   INQUIRE/KEY   ANSW/C/1/3  "Start anyway (y/n, default no) ?"
 
515
   IF ANSW(1:1) .EQ. "Y" THEN
 
516
           GOTO END
 
517
   ELSE
 
518
           WRITE/DESCR  middummc.tbl  TBLCONTR/I/4/1   0
 
519
           WRITE/DESCR  middummr.tbl  TBLCONTR/I/4/1   0
 
520
           GOTO  METH2
 
521
   ENDIF
 
522
ENDIF
 
523
 
 
524
END:
 
525
RETURN
 
526
 
 
527
! ************************ Start with a guess *************************
 
528
!
 
529
! The dispersion relation m*lambda = f(x,m) and order numbers
 
530
! m = f(x,y) have been stored (REGR, RORD) as polynomial regression
 
531
! coefficients. The different steps necessary to generate the
 
532
! solution for a new observation set include:
 
533
!
 
534
! - Determine the absolute order number of the lines in the table
 
535
!   using the solution RORD.
 
536
! - Determine the shift (due to flexures, etc...) by cross-correlation
 
537
!   of the X positions in the two tables using common orders.
 
538
! - Apply the X shift to the positions and generate an initial dispersion
 
539
!   relation.
 
540
!
 
541
ENTRY GUESS
 
542
!
 
543
DEFINE/LOCAL   ORDMAX/I/1/1   0    ? +lower 
 
544
DEFINE/LOCAL   SHIFLOG/I/1/1  0    ? +lower 
 
545
! Indicates if a shift value is provided in GUESS
 
546
DEFINE/LOCAL   SHIFT/R/1/1    0.   ? +lower 
 
547
DEFINE/LOCAL   GUETAB/C/1/60  XXX  ? +lower  ! Guess LINE table
 
548
DEFINE/LOCAL   MIDORD/I/1/2 0,0    ? +lower 
 
549
!
 
550
SHIFLOG = M$INDEX(GUESS,",")
 
551
IF SHIFLOG .GT. 0  THEN
 
552
   SHIFLOG = SHIFLOG - 1
 
553
   WRITE/KEY GUETAB  {GUESS(1:{SHIFLOG})}LINE.tbl
 
554
   SHIFLOG = SHIFLOG + 2
 
555
   SHIFT   = {GUESS({SHIFLOG}:)}
 
556
ELSE
 
557
   WRITE/KEY GUETAB  {GUESS}LINE.tbl
 
558
ENDIF
 
559
!
 
560
! Check existence of the guess table and copies the 
 
561
! average pixel value.
 
562
VERIFY/ECHELLE {GUETAB}  TAB
 
563
COPY/DK     {GUETAB} PIXEL/R/1/1    PIXEL/R/1/1  
 
564
IF {LWCL(1)} .GT. 0. THEN
 
565
        PIXEL(2) = LWCL(1) * PIXEL(1)
 
566
ELSE
 
567
        PIXEL(2) = M$ABS(LWCL(1))
 
568
ENDIF
 
569
!
 
570
! Determines absolute order numbers in the line.tbl
 
571
! and updates the ORDPOS keyword.
 
572
!
 
573
define/local doris/I/1/1 0
 
574
doris = M$EXISTD(GUETAB,"RORDI")
 
575
 
 
576
if doris .le. 0 then
 
577
   write/out "Warning:"
 
578
   write/out "Guess table {GUETAB} is of old format"
 
579
   write/out "Updating descriptors..."
 
580
   regress/poly {GUETAB} :ORDER :X,:YNEW {DEGREE},{DEGREE} KEYLONG {SESSOUTV}
 
581
   save/regr    {GUETAB} RORD  KEYLONG
 
582
endif
 
583
 
 
584
COPY/DD     {GUETAB} RORDI  {LINTAB}
 
585
COPY/DD     {GUETAB} RORDC  {LINTAB}
 
586
COPY/DD     {GUETAB} RORDR  {LINTAB}
 
587
COPY/DD     {GUETAB} RORDD  {LINTAB}
 
588
COMPUTE/REGR  {LINTAB} :AUX = RORD
 
589
COMPUTE/TABLE {LINTAB} :ORDER = INT(:AUX)
 
590
STATIST/TABLE {LINTAB} :ORDER {SESSOUTV}
 
591
ORDPOS(1) = outputr(2) ! Max
 
592
ordpos(2) = outputr(1) ! Min
 
593
echord(2) = ordpos(2)
 
594
echord(3) = ordpos(1)
 
595
!
 
596
! Determines the range of common absolute orders.
 
597
! Take the min of the two maximum absolute order numbers
 
598
! Take the max of the two minimum absolute order numbers
 
599
!
 
600
define/local guesord/I/1/2 0,0
 
601
define/local orange/I/1/2  0,0
 
602
COPY/DK  {GUETAB} ORDER/I/1/2  guesord
 
603
orange(1) = ordpos(1)
 
604
if guesord(1) .lt. orange(1)  orange(1) = guesord(1)
 
605
orange(2) = ordpos(2)
 
606
if guesord(2) .gt. orange(2)  orange(2) = guesord(2)
 
607
if orange(1) .lt. orange(2) then
 
608
   write/out "Guess Wavelength Calibration:"
 
609
   write/out "   No common absolute order number with session {GUESS}"
 
610
   write/out "   Cross-correlation is not peformed"
 
611
   shift = 0.
 
612
   goto gfit
 
613
endif
 
614
!
 
615
! Estimates the shift by cross-correlation or read value provided in 
 
616
! keyword GUESS. The estimate is performed on five orders selected in the 
 
617
! middle of the spectrum.
 
618
!
 
619
IF SHIFLOG .LE. 0 THEN
 
620
   MIDORD(1) = orange(2)
 
621
   MIDORD(2) = orange(1)
 
622
    WRITE/OUT "Estimate shift by cross-correlation"
 
623
    CORREL/LINE  {GUETAB}  {LINTAB}  ?  0.,0.25,100.,0.5 X,ORDER,+ {MIDORD(1)},{MIDORD(2)}  {SESSOUTV}
 
624
    SHIFT = OUTPUTR(1)
 
625
ENDIF
 
626
WRITE/OUT    "Cross-correlation shift value : 'SHIFT'"
 
627
!
 
628
! Applies the shift to the guess table and computes new coefficients
 
629
!
 
630
gfit:
 
631
COPY/TABLE      {GUETAB}  &g
 
632
COMPUTE/TABLE   &g :X = :X + 'shift'
 
633
COMPUTE/TABLE   &g :AUX = NULL
 
634
COMPUTE/TABLE   &g :AUX = :ORDER * :IDENT
 
635
SELECT/TABLE &g  (:SELECT.GT.0.5)  {SESSOUTV}
 
636
REGRES/POLY  &g  :AUX  :X,:ORDER  {DEGREE},{DEGREE}  KEYLONG {SESSOUTV}
 
637
SAVE/REGR    {LINTAB} REGR  KEYLONG
 
638
!
 
639
OUTPUTR(5) = {{LINTAB},REGRR(5)}
 
640
NBR(1)     = {{LINTAB},REGRI(1)}
 
641
RMS(1)     = 2.*OUTPUTR(5)/(ORDPOS(1)+ORDPOS(2))
 
642
!
 
643
WRITE/DESCR    {LINTAB}  ORDER/I/1/2   {ORDPOS(1)},{ORDPOS(2)}
 
644
WRITE/DESCR    {LINTAB}  PIXEL/R/1/2   {PIXEL(1)},{PIXEL(2)}
 
645
!
 
646
RETURN
 
647
! ************************  Load table in overlay plane *****************
 
648
 
 
649
ENTRY OVERLAY
 
650
 
 
651
DEFINE/LOCAL   NBROW/I/1/1   0    ? +lower 
 
652
 
 
653
IF WLCVISU(1:1) .NE. "Y"  RETURN
 
654
 
 
655
SELECT/TABLE   &r (sequence.ne.0)  {SESSOUTV}
 
656
NBROW = OUTPUTI(1)
 
657
COMPUTE/TABLE  &r  :SEQU = sequence
 
658
WRITE/DESCR    middummr.tbl  TBLCONTR/I/4/1  {NBROW}
 
659
CLEAR/CHAN    over
 
660
SELECT/TABLE  &r  (sequence.eq.1) {SESSOUTV}
 
661
LOAD/TABLE    &r  :X  :Y  :SEQU 0   10  {SESSOUTV}
 
662
SELECT/TABLE  &r  (sequence.ne.1) {SESSOUTV}
 
663
LOAD/TABLE    &r  :X  :Y  :SEQU 1   10 {SESSOUTV}
 
664
 
 
665
RETURN
 
666
 
 
667
! ************************  Cursor interaction **************************
 
668
 
 
669
ENTRY CURSOR
 
670
 
 
671
DEFINE/PARAMETER   P1    HAY   C    "Method" 
 
672
 
 
673
IF P1(1:1) .EQ. "R"  THEN
 
674
        IF OUTMODE(1:1) .EQ. "V" WRITE/OUT "Selection read from table middummr.tbl"
 
675
        WRITE/OUT "Previous method : {WLCMTD(1:8)}"
 
676
        RETURN
 
677
ENDIF
 
678
 
 
679
IF WLCVISU(1:1) .NE. "Y"  THEN
 
680
   @ creifnot 2 heat
 
681
   clear/chan over
 
682
   LOAD {WLC}
 
683
   WLCVISU = "YES"
 
684
ENDIF
 
685
 
 
686
! If the command CENTER/MOMENT is modified here, check consistency with
 
687
! command used near label RECLICK
 
688
 
 
689
WRITE/OUT    "For each line (for help, see CENTER/MOMENT with option CURSOR):"
 
690
WRITE/OUT    "       - select approximate position with cursor"
 
691
WRITE/OUT    "       - if necessary, optimize size of cursor rectangle"
 
692
WRITE/OUT    "       - enter position by pressing Enter on keyboard"
 
693
WRITE/OUT    "Exit by pressing Exit button"
 
694
 
 
695
CENTER/MOMENT  CURSOR  &c  
 
696
 
 
697
COMPUTE/TABLE  middummc   :XERR = ABS(:XEND - :XSTART)/2.
 
698
COMPUTE/TABLE  middummc   :YERR = ABS(:YEND - :YSTART)/2.
 
699
 
 
700
SELECT/TABLE   middummc   (ABS(:STATUS).LT.0.01)  {SESSOUTV}
 
701
 
 
702
COPY/TT        middummc   :XCEN     middummr   :X  {SESSOUTV}
 
703
COPY/TT        middummc   :YCEN     middummr   :Y  {SESSOUTV}
 
704
COPY/TT        middummc   :XERR     middummr   :DX {SESSOUTV}
 
705
COPY/TT        middummc   :YERR     middummr   :DY {SESSOUTV}
 
706
 
 
707
WRITE/DESCR    middummr.tbl    WLCMTD/C/1/8   {METHOD(1:5)}
 
708
 
 
709
RETURN
 
710
 
 
711
! ********************* Prepares association loop values *******************
 
712
 
 
713
ENTRY PREPLOOP
 
714
 
 
715
COMPUTE/TABLE  {LINTAB}   :IDENT = 0.
 
716
 
 
717
LAST(1) = {middummr.tbl,TBLCONTR(4)}
 
718
 
 
719
RETURN
 
720
 
 
721
! ******************** Associates positions ****************************
 
722
 
 
723
ENTRY ASSOC
 
724
 
 
725
STASS:
 
726
 
 
727
     xypos(1) = {middummr.tbl,:X,@{line}}
 
728
     xypos(2) = {middummr.tbl,:Y,@{line}}
 
729
     xymax(1) = {middummr.tbl,:DX,@{line}}
 
730
     xymax(2) = {middummr.tbl,:DY,@{line}}
 
731
 
 
732
     select/table  {LINTAB} ABS(:YNEW-'xypos(2)').LT.'xymax(2)'.and.ABS(:X-'xypos(1)').LT.'xymax(1)' {SESSOUTV}
 
733
     nbsel = outputi(1)
 
734
 
 
735
     IF nbsel .ne. 1 THEN
 
736
       WRITE/OUT "*** Sorry. Cannot identify feature {LINE}        ***"
 
737
       WRITE/OUT "*** Line either not found or blended             ***"
 
738
       WRITE/OUT "*** Please start over again with new line        ***"
 
739
       WRITE/OUT "*** All detected lines are now marked on display ***"
 
740
       SELECT/TABLE  {LINTAB}  ALL
 
741
       LOAD/TABLE    {LINTAB}  :X  :YNEW  {SESSOUTV}
 
742
       @s  neciden,reenter
 
743
       GOTO STASS
 
744
     ENDIF
 
745
 
 
746
     COPY/TT  {LINTAB}   :SEQU   middummt   #1 {SESSOUTV}
 
747
     SEQU(1)   =  {middummt.tbl,#1,@1}
 
748
     COPY/TT  {LINTAB}   :Y      middummt   #1 {SESSOUTV}
 
749
     ORDVAL(1)   =  {middummt.tbl,#1,@1}
 
750
     WRITE/TAB  {LINTAB} :SELECT @{SEQU}  2
 
751
 
 
752
     IF line .eq. 1  THEN
 
753
       INQUIRE/KEYW COUPL/I/1/1  "Enter absolute order number of first pointed line (square mark) : "
 
754
       ORDIFF(1)  =  M$NINT(ORDVAL(1)+COUPL(1))
 
755
       COMPUTE/TABLE  {LINTAB}   :ORDER = INT('ORDIFF'-:Y) !Nearest integer
 
756
       ORDPOS(1)  =   ORDIFF  -  NBORD
 
757
       ORDPOS(2)  =   ORDIFF  -  1
 
758
       WRITE/DESCR    {LINTAB}    ORDER/I/1/2    'ORDPOS(2)','ORDPOS(1)'
 
759
       ECHORD(3)  =   ORDPOS(2)   ! Non desirable duplication of info.
 
760
       ECHORD(2)  =   ORDPOS(1)
 
761
     ENDIF
 
762
 
 
763
RETURN
 
764
 
 
765
! ******************************************************************
 
766
 
 
767
ENTRY ITER
 
768
 
 
769
DEFINE/LOCAL ITER/I/1/1   0              ? +lower 
 
770
DEFINE/LOCAL NITER/I/1/1  {WLCNITER(2)}  ? +lower 
 
771
DEFINE/LOCAL NINDEP/I/1/1 0              ? +lower 
 
772
DEFINE/LOCAL PASS/I/1/1   0              ? +lower 
 
773
define/local rmsratio/R/1/1  0.
 
774
 
 
775
SET/FORMAT F12.5  I1
 
776
 
 
777
WRITE/OUT "Iterative Identification:"
 
778
 
 
779
ITERAT:
 
780
 
 
781
COMPUTE/REGRES  {LINTAB}  :WAVEC =  REGR
 
782
IF WLCMTD(1:1) .NE. "C"  THEN
 
783
   COMPUTE/TABLE   {LINTAB}  :WAVEC = :WAVEC/:ORDER
 
784
ENDIF
 
785
COMPUTE/TABLE   {LINTAB}  :IDENT = 0.
 
786
 
 
787
WRITE/KEY   IN_A/C/1/60     {LINTAB}
 
788
WRITE/KEY   IN_B/C/1/60     {LINCAT}
 
789
WRITE/KEY   INPUTR/R/1/2    {LWCL(2)}
 
790
WRITE/KEY   OUT_A/C/1/60    middummk.bdf
 
791
WRITE/KEY   INPUTI/I/1/1    400
 
792
WRITE/KEY   INPUTD/D/1/2    -5.,0.05
 
793
inputr(2) = 1.
 
794
 
 
795
SET/MIDAS OUTPUT=NO
 
796
WRITE/KEY   PROGSTAT/I/1/1   0
 
797
RUN STD_EXE:spematch  
 
798
IF PROGSTAT(1) .NE. 0  GOTO EXIT ! Problems in matching loop
 
799
SET/MIDAS OUTPUT=YES
 
800
 
 
801
!KB SELECT/TAB       {LINTAB} (:IDENT.NE.0.0)  {SESSOUTV}
 
802
SELECT/TAB       {LINTAB} (ABS(:IDENT).GT.0.00001)  {SESSOUTV}
 
803
 
 
804
NBR(2) = OUTPUTI(1)
 
805
WRITE/OUT  "*** Number of identifications: {OUTPUTI(1)}"
 
806
 
 
807
IF WLCVISU(1:1) .EQ. "Y" THEN
 
808
   CLEAR/CHAN       OVER
 
809
   LOAD/TAB         {LINTAB}  :X  :YNEW  {SESSOUTV}
 
810
ENDIF
 
811
 
 
812
STAT/TABLE       {LINTAB}  :RESIDUAL  {SESSOUTV}
 
813
RMS(2)   = OUTPUTR(4)
 
814
PIXEL(3) = M$ABS(RMS(2)/PIXEL(1))
 
815
 
 
816
! Tests whether one stops the iteration
 
817
IF ITER .EQ. NITER  GOTO ENDIT
 
818
ITER = ITER + 1
 
819
IF ITER .GT. 1 THEN
 
820
  rmsratio = M$ABS(1. - rms(2)/rms(1))
 
821
  IF NBR(2) .LT. NBR(1) .OR. RMSRATIO .LT. 0.01 THEN
 
822
     IF ITER .GT. WLCNITER(1) GOTO  ENDIT
 
823
  ENDIF
 
824
ENDIF
 
825
NBR(1) = NBR(2)
 
826
RMS(1) = RMS(2)
 
827
 
 
828
WRITE/OUT  "*** STD. DEV. {RMS(2)} wavelength units, i.e. {PIXEL(3)} pixels  ***"
 
829
 
 
830
IF PIXEL(3) .GT. LWCL(3)  THEN
 
831
   WRITE/OUT  "Sorry, wrong identifications..."
 
832
   RETURN/EXIT
 
833
ENDIF
 
834
 
 
835
COMPUTE/TABLE   {LINTAB}   :AUX = 0.
 
836
COMPUTE/TABLE   {LINTAB}   :AUX = :ORDER * :IDENT
 
837
!KB SELECT/TABLE    {LINTAB}   (:IDENT.NE.0.0) {SESSOUTV}
 
838
SELECT/TAB       {LINTAB} (ABS(:IDENT).GT.0.00001)  {SESSOUTV}
 
839
 
 
840
 
 
841
IF NBR(2) .GT. NBRMIN THEN
 
842
     WRITE/KEY  WLCMTD  TWO-D
 
843
ENDIF
 
844
 
 
845
IF WLCMTD(1:1) .EQ. "P"  REGRES/POLY  {LINTAB} :AUX :XROT      {DEGREE} KEYLONG {SESSOUTV}
 
846
IF WLCMTD(1:1) .EQ. "A"  REGRES/POLY  {LINTAB} :AUX :XROT      {DEGREE} KEYLONG {SESSOUTV}
 
847
 
 
848
IF WLCMTD(1:1) .EQ. "L"  REGRES/POLY  {LINTAB} :AUX :X,:ORDER  {DEGREE},{DEGREE} KEYLONG {SESSOUTV}
 
849
IF WLCMTD(1:1) .EQ. "G"  REGRES/POLY  {LINTAB} :AUX :X,:ORDER  {DEGREE},{DEGREE} KEYLONG {SESSOUTV}
 
850
IF WLCMTD(1:1) .EQ. "T"  REGRES/POLY  {LINTAB} :AUX :X,:ORDER  {DEGREE},{DEGREE} KEYLONG {SESSOUTV}
 
851
 
 
852
IF WLCMTD(1:1) .EQ. "C"  REGRES/POLY  {LINTAB} :IDENT  :YNEW   {DEGREE} KEYLONG {SESSOUTV}
 
853
 
 
854
SAVE/REGRES       {LINTAB}   REGR  KEYLONG
 
855
 
 
856
GOTO ITERAT
 
857
 
 
858
ENDIT:
 
859
 
 
860
SET/FORMAT
 
861
RETURN
 
862
 
 
863
!*****************************************************************************
 
864
 
 
865
ENTRY ANALYS
 
866
 
 
867
   SELECT/TABLE  {LINTAB}  (:SELECT.GT.1) {SESSOUTV}
 
868
   WRITE/KEY  IN_A  {LINTAB}
 
869
   WRITE/KEY  IN_B  {WLCMTD}
 
870
   inputr(2) = 1.
 
871
   RUN STD_EXE:necwlcopt
 
872
   IF OUTPUTI(1) .NE. 0  RETURN/EXIT
 
873
 
 
874
RETURN
 
875
 
 
876
!************************************************************************
 
877
 
 
878
ENTRY REENTER
 
879
 
 
880
IF WLCMTD(1:1) .EQ. "P" THEN
 
881
  WRITE/OUT "Please identify the pair again, in the same order"
 
882
  WRITE/OUT "as you did it the first time."
 
883
ELSE
 
884
  WRITE/OUT "Please identify line {line} again."
 
885
ENDIF
 
886
 
 
887
CENTER/MOMENT CURSOR &c
 
888
 
 
889
COMPUTE/TABLE  middummc   :XERR = ABS(:XEND - :XSTART)/2.
 
890
COMPUTE/TABLE  middummc   :YERR = ABS(:YEND - :YSTART)/2.
 
891
 
 
892
SELECT/TABLE   middummc   (ABS(:STATUS).LT.0.01)  {SESSOUTV}
 
893
 
 
894
COPY/TAB       &c   &d
 
895
 
 
896
DEFINE/LOCAL   LINPAIR/I/1/1  0   ? +lower 
 
897
IF line .le. 2 THEN
 
898
          linpair = 1
 
899
ELSE
 
900
          linpair = 3
 
901
ENDIF
 
902
 
 
903
IF WLCMTD(1:1) .EQ. "P" THEN
 
904
   WRITE/TABLE &r :X  @{linpair}  {&d,:XCEN,@1}   
 
905
   WRITE/TABLE &r :Y  @{linpair}  {&d,:YCEN,@1}   
 
906
   WRITE/TABLE &r :DX @{linpair}  {&d,:XERR,@1}   
 
907
   WRITE/TABLE &r :DY @{linpair}  {&d,:YERR,@1}   
 
908
   linpair = linpair + 1
 
909
   WRITE/TABLE &r :X  @{linpair}  {&d,:XCEN,@2}   
 
910
   WRITE/TABLE &r :Y  @{linpair}  {&d,:YCEN,@2}   
 
911
   WRITE/TABLE &r :DX @{linpair}  {&d,:XERR,@2}   
 
912
   WRITE/TABLE &r :DY @{linpair}  {&d,:YERR,@2}   
 
913
ELSE
 
914
   WRITE/TABLE &r :X  @{line}  {&d,:XCEN,@1}   
 
915
   WRITE/TABLE &r :Y  @{line}  {&d,:YCEN,@1}   
 
916
   WRITE/TABLE &r :DX @{line}  {&d,:XERR,@1}   
 
917
   WRITE/TABLE &r :DY @{line}  {&d,:YERR,@1}   
 
918
ENDIF
 
919
 
 
920
RETURN
 
921
 
 
922
 
 
923
!***********************************************************************
 
924
 
 
925
ENTRY HOUGH
 
926
 
 
927
define/local    factor/D/1/1  0.
 
928
define/local    convert/I/1/1 0
 
929
define/local    ymean/R/1/1   0.
 
930
define/local    rms/R/1/1     0.
 
931
define/local    slop/D/1/1    0.
 
932
 
 
933
convert = relord + absord
 
934
 
 
935
echord(6) = convert
 
936
echord(2) = convert - echord(1) 
 
937
echord(3) = convert - 1
 
938
 
 
939
ordpos(1) = echord(3)
 
940
ordpos(2) = echord(2)
 
941
 
 
942
copy/kd ordpos/I/1/2 {lintab} ORDER/I/1/2
 
943
 
 
944
compute/table {LINTAB}  :ORDER = {CONVERT} - :Y
 
945
 
 
946
SELECT/TABLE {LINTAB} :Y.EQ.{relord} {SESSOUTV}
 
947
 
 
948
stat/table   {LINTAB} :YNEW  {SESSOUTV}
 
949
ymean = outputr(3)
 
950
 
 
951
dispersion/hough ? ? {IMSIZE(1)},{{WLC},START(1)},{{WLC},STEP(1)}  {LINTAB}  {SESSOUTV}
 
952
 
 
953
COMPUTE/TABLE   {LINTAB}  :AUX  = :WAVEC*:ORDER
 
954
SELECT/TABLE    {LINTAB}  :Y.EQ.{relord}  {SESSOUTV}
 
955
REGRESS/POLY    {LINTAB}  :AUX  :X  1     {SESSOUTV}
 
956
 
 
957
slop = outputd(2)
 
958
rms = lwcl(4)*slop
 
959
rms = rms*rms
 
960
 
 
961
! First estimate on the error on the angle = 1. degree
 
962
factor = slop*M$TAN(1)
 
963
 
 
964
SAVE/REGR       {LINTAB}  REGR
 
965
 
 
966
COMPUTE/TABLE   {LINTAB}  :ERROR = (:YNEW-{YMEAN})*{FACTOR}
 
967
COMPUTE/TABLE   {LINTAB}  :ERROR = SQRT(:ERROR*:ERROR + {RMS})
 
968
COMPUTE/TABLE   {LINTAB}  :ERROR = :ERROR/:ORDER
 
969
 
 
970
stat/tab {lintab} :ERROR  {SESSOUTV}
 
971
write/out  "Error range: {outputr(1)} to {outputr(2)} wav. units"
 
972
 
 
973
SELECT/TABLE {LINTAB} ALL
 
974
pixel(1) = avdisp
 
975
 
 
976
RETURN
 
977
 
 
978
! ******************************************************************
 
979
 
 
980
ENTRY HTITER
 
981
 
 
982
DEFINE/LOCAL ITER/I/1/1   0              ? +lower 
 
983
DEFINE/LOCAL NITER/I/1/1  {WLCNITER(2)}  ? +lower 
 
984
DEFINE/LOCAL NINDEP/I/1/1 0              ? +lower 
 
985
DEFINE/LOCAL PASS/I/1/1   0              ? +lower 
 
986
define/local rmsratio/R/1/1  0.
 
987
 
 
988
WRITE/OUT ".. Iterative Identification:"
 
989
 
 
990
IF WLCVISU(1:1) .EQ. "Y" THEN
 
991
   graph/spec
 
992
   set/graph
 
993
   set/graph pmode=1
 
994
ENDIF
 
995
 
 
996
HTITER:
 
997
 
 
998
COMPUTE/REGRES  {LINTAB}  :WAVEC =  REGR
 
999
COMPUTE/TABLE   {LINTAB}  :WAVEC = :WAVEC/:ORDER
 
1000
 
 
1001
SET/FORMAT F12.5  I1
 
1002
 
 
1003
COMPUTE/TABLE   {LINTAB}  :IDENT  = 0.
 
1004
COMPUTE/TABLE   {LINTAB}  :SELECT = 0.
 
1005
 
 
1006
WRITE/KEY   IN_A/C/1/60     {LINTAB}
 
1007
WRITE/KEY   IN_B/C/1/60     {LINCAT}
 
1008
WRITE/KEY   INPUTR/R/1/2    {LWCL(2)}
 
1009
WRITE/KEY   OUT_A/C/1/60    middummk.bdf
 
1010
WRITE/KEY   INPUTI/I/1/1    400
 
1011
WRITE/KEY   INPUTD/D/1/2    -5.,0.05
 
1012
inputr(2) = 1.
 
1013
 
 
1014
SET/MIDAS OUTPUT=NO
 
1015
WRITE/KEY   PROGSTAT/I/1/1   0
 
1016
RUN STD_EXE:spematch  
 
1017
IF PROGSTAT(1) .NE. 0  GOTO EXIT ! Problems in matching loop
 
1018
SET/MIDAS OUTPUT=YES
 
1019
 
 
1020
STAT/TABLE {LINTAB} :RESIDUAL
 
1021
RMS(2)   = OUTPUTR(4)*2./(echord(2)+echord(3))
 
1022
PIXEL(3) = RMS(2)/PIXEL(1)
 
1023
WRITE/OUT  "*** Matching: STD. DEV. {RMS(2)} wavelength units, i.e. {PIXEL(3)} pixels  ***"
 
1024
 
 
1025
IF WLCVISU(1:1) .EQ. "Y" THEN
 
1026
  PLOT &k
 
1027
  PLOT/HISTO {LINTAB} :RESIDUAL 
 
1028
  select/table  {LINTAB} :IDENT.GT.0.0
 
1029
  PLOT/TABLE {LINTAB} :IDENT :RESIDUAL
 
1030
ENDIF
 
1031
 
 
1032
!KB SELECT/TAB       {LINTAB} (:IDENT.NE.0.0)  {SESSOUTV}
 
1033
SELECT/TAB       {LINTAB} (ABS(:IDENT).GT.0.00001)  {SESSOUTV}
 
1034
 
 
1035
NBR(2) = OUTPUTI(1)
 
1036
WRITE/OUT  "*** Iteration: {ITER}. Number of identifications: {OUTPUTI(1)}"
 
1037
 
 
1038
IF WLCVISU(1:1) .EQ. "Y" THEN
 
1039
   !KB SELECT/TAB       {LINTAB} (:IDENT.NE.0.0)  {SESSOUTV}
 
1040
   SELECT/TAB       {LINTAB} (ABS(:IDENT).GT.0.00001)  {SESSOUTV}
 
1041
 
 
1042
   CLEAR/CHAN       OVER
 
1043
   LOAD/TAB         {LINTAB}  :X  :YNEW  {SESSOUTV}
 
1044
ENDIF
 
1045
 
 
1046
COMPUTE/TABLE   {LINTAB}   :AUX = :ORDER * :IDENT
 
1047
SELECT/TAB      {LINTAB} (:IDENT.GT.0.0)  {SESSOUTV}
 
1048
 
 
1049
IF NBR(2) .GT. NBRMIN THEN
 
1050
     WRITE/KEY  WLCMTD  TWO-D
 
1051
ENDIF
 
1052
 
 
1053
 
 
1054
! Robust Regression PROGRESS
 
1055
compute/table   {LINTAB} :X2  = :X*:X
 
1056
compute/table   {LINTAB} :X3  = :X*:X*:X
 
1057
compute/table   {LINTAB} :O2  = :ORDER*:ORDER
 
1058
compute/table   {LINTAB} :O3  = :O2*:ORDER
 
1059
compute/table   {LINTAB} :O4  = :O2*:O2
 
1060
compute/table   {LINTAB} :O3X = :O3*:X
 
1061
compute/table   {LINTAB} :O3X2= :O3*:X2
 
1062
compute/table   {LINTAB} :OX  = :ORDER*:X
 
1063
compute/table   {LINTAB} :O2X  = :O2*:X
 
1064
compute/table   {LINTAB} :OX2  = :ORDER*:X2
 
1065
compute/table   {LINTAB} :O2X2 = :O2*:X2
 
1066
!compute/table   {LINTAB} :PIXFRAC = :X - INT(:X)
 
1067
 
 
1068
select/table    {LINTAB} :IDENT.GT.0.0 {SESSOUTV}
 
1069
regress/robust  {LINTAB} :AUX :ORDER,:X,:X2,:OX,:O2,:O2X,:O3,:O3X,:O4,:X3 ? :AUX  {SESSOUTV}
 
1070
 
 
1071
!Computes :WAVEC for all cases
 
1072
!set/format F15.5
 
1073
!COMPUTE/TABLE {LINTAB} :WAVEC = {OUTPUTR(1)}*:ORDER+ {OUTPUTR(2)}*:X+ {OUTPUTR(3)}*:X2+ {OUTPUTR(4)}*:X3
 
1074
!COMPUTE/TABLE {LINTAB} :WAVEC = :WAVEC + {OUTPUTR(5)}*:OX + {OUTPUTR(6)}*:O2 + {OUTPUTR(7)}*:O2X
 
1075
!COMPUTE/TABLE {LINTAB} :WAVEC = :WAVEC + {OUTPUTR(8)}*:OX2 + {OUTPUTR(9)}*:O2X2 + {OUTPUTR(10)}
 
1076
 
 
1077
! Determines median of residuals
 
1078
 
 
1079
select/table  {lintab} :IDENT.GT.0.0
 
1080
 
 
1081
IF WLCVISU(1:1) .EQ. "Y" THEN
 
1082
  set/graph colour=1 pmode=1
 
1083
  PLOT/TABLE {LINTAB}  :IDENT :RESIDUAL
 
1084
ENDIF
 
1085
 
 
1086
copy/table    {LINTAB} &l
 
1087
compute/table &l :R2 = ABS(:RESIDUAL)
 
1088
sort/table    &l :R2
 
1089
inputi(1)  = {middumml.tbl,TBLCONTR(4)}/2
 
1090
outputr(4) = {middumml.tbl,:R2,@{inputi(1)}}
 
1091
RMS(2)     = OUTPUTR(4)*2.5 
 
1092
 
 
1093
SELECT/TABLE {LINTAB} :IDENT.GT.0.0.AND.ABS(:RESIDUAL).LT.{RMS(2)} {SESSOUTV}
 
1094
 
 
1095
IF WLCVISU(1:1) .EQ. "Y" THEN
 
1096
  select/table {lintab} .NOT.SELECT
 
1097
  set/graph colour=4
 
1098
  OVER/TABLE {lintab} :IDENT :RESIDUAL
 
1099
  select/table {lintab} .NOT.SELECT
 
1100
  set/graph colour=0
 
1101
ENDIF
 
1102
 
 
1103
 
 
1104
! Measures filtered dispersion 
 
1105
STAT/TABLE   {LINTAB} :RESIDUAL {SESSOUTV}
 
1106
RMS(2)   = OUTPUTR(4)*2./(echord(2)+echord(3))
 
1107
COMPUTE/TABLE   {LINTAB}    :ERROR = {OUTPUTR(4)}
 
1108
COMPUTE/TABLE   {LINTAB}    :ERROR = :ERROR/:ORDER
 
1109
PIXEL(3) = RMS(2)/PIXEL(1)
 
1110
WRITE/OUT  "*** Rejection: STD. DEV. {RMS(2)} wav. units, i.e. {PIXEL(3)} pixels  ***"
 
1111
 
 
1112
!Computes wavelength by polynomial interpolation
 
1113
!IF WLCMTD(1:1) .EQ. "P"  REGRES/POLY  {LINTAB} :AUX :XROT      {DEGREE} KEYLONG {SESSOUTV}
 
1114
!IF WLCMTD(1:1) .EQ. "A"  REGRES/POLY  {LINTAB} :AUX :XROT      {DEGREE} KEYLONG {SESSOUTV}
 
1115
!IF WLCMTD(1:1) .EQ. "L"  REGRES/POLY  {LINTAB} :AUX :X,:ORDER  {DEGREE},{DEGREE} KEYLONG {SESSOUTV}
 
1116
!IF WLCMTD(1:1) .EQ. "G"  REGRES/POLY  {LINTAB} :AUX :X,:ORDER  {DEGREE},{DEGREE} KEYLONG {SESSOUTV}
 
1117
!IF WLCMTD(1:1) .EQ. "T"  REGRES/POLY  {LINTAB} :AUX :X,:ORDER  {DEGREE},{DEGREE} KEYLONG {SESSOUTV}
 
1118
 
 
1119
REGRES/POLY  {LINTAB} :AUX :X,:ORDER  {DEGREE},{DEGREE} KEYLONG {SESSOUTV}
 
1120
SAVE/REGR    {LINTAB} REGR  KEYLONG {SESSOUTV} 
 
1121
 
 
1122
! Tests whether one stops the iteration
 
1123
IF ITER .GT. NITER  THEN
 
1124
   WRITE/OUT "Reached maximal number of iterations WLCNITER(2)={WLCNITER(2)}"
 
1125
   GOTO HTENDIT
 
1126
ENDIF
 
1127
 
 
1128
ITER = ITER + 1
 
1129
IF ITER .GT. 1 THEN
 
1130
  rmsratio = M$ABS(1. - rms(2)/rms(1))
 
1131
  IF NBR(2) .LT. NBR(1) .OR. RMSRATIO .LT. 0.01 THEN
 
1132
       IF ITER .GT. WLCNITER(1) GOTO  HTENDIT
 
1133
  ENDIF
 
1134
ENDIF
 
1135
NBR(1) = NBR(2)
 
1136
RMS(1) = RMS(2)
 
1137
 
 
1138
 
 
1139
IF PIXEL(3) .GT. LWCL(3)  THEN
 
1140
   WRITE/OUT  "Sorry, wrong identifications..."
 
1141
   RETURN/EXIT
 
1142
ENDIF
 
1143
 
 
1144
GOTO HTITER
 
1145
 
 
1146
HTENDIT:
 
1147
 
 
1148
SET/FORMAT
 
1149
 
 
1150
RETURN
 
1151
 
 
1152
 
 
1153
!***************************************************************************
 
1154
 
 
1155
ENTRY VERIF
 
1156
 
 
1157
DEFINE/LOCAL RATIO/R/1/1  0. ? +lower 
 
1158
DEFINE/LOCAL BRIGHT/I/1/1 0  ? +lower 
 
1159
DEFINE/LOCAL SELECT/I/1/1 0  ? +lower 
 
1160
DEFINE/LOCAL MINSEL/I/1/1 0  ? +lower 
 
1161
 
 
1162
SELECT/TABLE {LINTAB}  ALL
 
1163
COPY/TI      {LINTAB}  &h
 
1164
STATIST/IMA  &h  [<,@3:>,@3] OPTION=HN
 
1165
 
 
1166
SELECT/TABLE {LINTAB} :PEAK.GT.{OUTPUTR(8)} {SESSOUTV}  ! OUTPUTR(8) = Median
 
1167
BRIGHT = OUTPUTI(1)
 
1168
 
 
1169
SELECT/TABLE {LINTAB} SELECT.AND.:SELECT.GE.1 {SESSOUTV}
 
1170
SELECT = OUTPUTI(1)
 
1171
 
 
1172
SELECT/TABLE {LINTAB}  ALL
 
1173
MINSEL = DEGREE + 2
 
1174
RATIO  = 100.*SELECT/BRIGHT
 
1175
 
 
1176
SET/FORMAT  I1  F3.0
 
1177
 
 
1178
WRITE/OUT "**************************  Verification  ***********************************"
 
1179
WRITE/OUT "     "
 
1180
WRITE/OUT "1) Minimum number of selections per order : {MINSEL}"
 
1181
WRITE/OUT "   If the number of selections in any order (column NO.LINES above)"
 
1182
WRITE/OUT "   is less or equal than the minimum, this order should be checked."
 
1183
WRITE/OUT "        "
 
1184
WRITE/OUT "2) Percentage of identifications among the half brighter lines : {RATIO} %"
 
1185
WRITE/OUT "   This percentage must be as high as possible (above 50%). Low values"
 
1186
WRITE/OUT "   indicate an uncertain calibration."
 
1187
WRITE/OUT "        "
 
1188
 
 
1189
SET/FORMAT
 
1190
RETURN
 
1191
 
 
1192
 
 
1193
!**********************************************************************
 
1194
 
 
1195
ENTRY MAPRES
 
1196
 
 
1197
compute/table {LINTAB} :RES=:RESIDUAL * :RESIDUAL
 
1198
select/table {LINTAB} :IDENT.GT.0.0
 
1199
convert/table &l = line.tbl :X,:YNEW :RES {WLC} POLY 1,1
 
1200
 
 
1201
select/table {LINTAB}  all
 
1202
 
 
1203
COMPUTE/TABLE  {LINTAB}  :XSTA = :X - 1
 
1204
COMPUTE/TABLE  {LINTAB}  :XEND = :X + 1
 
1205
COMPUTE/TABLE  {LINTAB}  :YSTA = :YNEW - 1
 
1206
COMPUTE/TABLE  {LINTAB}  :YEND = :YNEW + 1
 
1207
COMPUTE/TABLE  {LINTAB}  :YBKG = :YNEW
 
1208
 
 
1209
WRITE/KEY IN_A/C/1/60  {LINTAB}
 
1210
WRITE/KEY IN_B/C/1/60  middumml.bdf
 
1211
WRITE/KEY OUT_A/C/1/60 middummi.bdf
 
1212
 
 
1213
WRITE/KEY INPUTI/I/1/2 1,-10
 
1214
 
 
1215
RUN STD_EXE:NECBACK
 
1216
 
 
1217
compute/table {LINTAB} :ERROR = 2.*SQRT(:BKG)
 
1218
 
 
1219
dele/col {lintab} :XSTA
 
1220
dele/col {lintab} :YSTA
 
1221
dele/col {lintab} :XEND
 
1222
dele/col {lintab} :YEND
 
1223
dele/col {lintab} :BKG
 
1224
dele/col {lintab} :RES
 
1225
 
 
1226
RETURN
 
1227