1
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++
2
!.COPYRIGHT (C) 1991-2011 European Southern Observatory
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
12
!-------------------------------------------------------
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}
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 "
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"
39
VERIFY/ECHELLE {LINTAB} TAB
40
VERIFY/SPEC {LINCAT} MID_ARC LINCAT TAB
42
IF WLCVISU(1:1) .EQ. "Y" THEN
47
WRITE/KEY IN_A {LINTAB}
48
WRITE/KEY IN_B {LINCAT}
49
INPUTR(1) = -10 ! Verification mode
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
59
IF INPUTR(1) .LT. 0 RETURN/EXIT
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
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
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
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
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
106
WRITE/OUT "*** Wavelength calibration ***"
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}"
128
IF LWCL(2) .LE. 0. ERROR/ECHELLE IDENTIFY/ECHELLE LWCL(2)
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.
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"
142
IF WLCMTD(1:1) .EQ. "O" THEN
143
COPY/DK {LINTAB} PIXEL PIXEL
147
WRITE/DESCR {LINTAB} HISTORY_UPDA/I/1/1 0 ! No update of history
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}
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}
161
CREATE/COLUMN &t :COL
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.
169
IF WLCMTD(1:1) .EQ. "R" THEN
170
COPY/DK middummr.tbl WLCMTD WLCMTD
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
183
IF WLCMTD(1:1) .EQ. "P" @s neciden,pair
185
IF WLCMTD(1:1) .EQ. "A" @s neciden,angle
186
IF WLCMTD(1:1) .EQ. "T" @s neciden,angle
188
IF WLCMTD(1:1) .EQ. "G" @s neciden,guess
189
IF WLCMTD(1:1) .EQ. "L" @s neciden,hough
191
COPY/KD PIXEL {LINTAB} PIXEL
192
COPY/KD WLCMTD middummr.tbl WLCMTD
194
! Branch to the main iteration loop
196
IF WLCREG(1:1) .EQ. "S" @s neciden,iter
197
IF WLCREG(1:1) .EQ. "R" @s neciden,htiter
200
pixel(1) = M$ABS(pixel(1))
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).
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
214
PIXEL(2) = 3.*KEYLONGR(5)*2./(echord(2)+echord(3))
217
PIXEL(2) = TOL*PIXEL(1)
219
PIXEL(2) = M$ABS(TOL)
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}
228
Write/Out "Final Selection: Threshold = {pixel(2)} wav. units, Nb = {outputi(1)}"
229
IF WLCVISU(1:1) .EQ. "Y" THEN
231
LOAD/TABLE {LINTAB} :X :YNEW
234
! Save the final dispersion relation and order numbers for later
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
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)}
245
INPUTI(2) = echord(1)
246
INPUTI(3) = imsize(1)
247
INPUTI(4) = {WLCOPT(1:1)}
249
INPUTR(1) = PIXEL(2) ! Tolerance in wavelength units
251
INPUTD(1) = {{WLC},START(1)}
252
INPUTD(2) = {{WLC},STEP(1)}
256
WLRANGE(1) = OUTPUTR(1)
257
WLRANGE(2) = OUTPUTR(2)
261
! Estimate value of K for optional IUE-like RIPPLE correction
263
STAT/TAB {LINTAB} :AUX {SESSOUTV}
264
write/key SAMPLE {{LINTAB},PIXEL(1)}
266
echord(6) = echord(3) + 1
267
write/out "Set parameters SAMPLE={SAMPLE} and RIPK={RIPK}"
271
if wlcvisu(1:1) .eq. "Y" PLOT/RESIDUAL
275
!*****************************************************************************
281
IF METHOD(1:1) .NE. "R" THEN
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."
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"
294
@s neciden,cursor {METHOD}
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
303
! *** Find selections in table LINE
314
COUPL(1) = {{LINTAB},:ORDER,@{SEQU}}
315
INQUIRE/KEY WAVE/D/1/1 "Sequence no. {LINE}, Order no. {COUPL}. Enter wavelength : "
318
COPY/KT WAVE/D/1/1 {LINTAB} :IDENT @{SEQU}
322
! computes rotation angle and peforms rotation
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})
328
ALPHA(1) = M$ATAN (ALPHA(1))
329
ALPHA(2) = M$ATAN (ALPHA(2))
332
! @s echidenew,analys
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."
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
353
COMPUTE/TABLE {LINTAB} :XROT = :X*{CCDBIN(1)}*COS('alpha(3)')+:YNEW*{CCDBIN(2)}*SIN('alpha(3)')
356
COMPUTE/TABLE {LINTAB} :AUX = :ORDER * :IDENT
357
SELECT/TABLE {LINTAB} (:IDENT.GT.0.0) {SESSOUTV}
359
REGRES/POLY {LINTAB} :AUX :XROT 2 {SESSOUTV}
360
SAVE/REGRES {LINTAB} REGR
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)
366
IF {LWCL(1)} .GT. 0. THEN
367
PIXEL(2) = LWCL(1)*PIXEL(1)
369
PIXEL(2) = M$ABS(LWCL(1))
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"
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."
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
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
397
WRITE/OUT "Inaccuracy is too large. Aborted..."
403
!Prepares new table with analytical rotation angle
404
COMPUTE/TABLE {LINTAB} :XROT = :X*{CCDBIN(1)}*COS('alpmax')+:YNEW*{CCDBIN(2)}*SIN('alpmax')
406
COMPUTE/TABLE {LINTAB} :AUX = :ORDER * :IDENT
407
SELECT/TABLE {LINTAB} (:IDENT.GT.0.0)
409
REGRES/POLY {LINTAB} :AUX :XROT 2
410
SAVE/REGRES {LINTAB} REGR
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))
416
IF {LWCL(1)} .GT. 0. THEN
417
PIXEL(2) = {LWCL(1)} * PIXEL(1)
419
PIXEL(2) = M$ABS({LWCL(1)})
424
! ***************** Start with 4 lines and play with rotation ************
428
DEFINE/LOCAL DEGSTART/I/1/1 2 ? +lower
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
438
IF METHOD(1:1) .NE. "R" THEN
439
WRITE/OUT "Point at least {NBLMIN} lines of which you know the wavelength"
442
@s neciden,cursor {METHOD}
445
outputi(1) = outputi(1) + 1
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
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}
466
!KB SELECT/TABLE {LINTAB} (:IDENT.NE.0.0) {SESSOUTV}
467
SELECT/TAB {LINTAB} (ABS(:IDENT).GT.0.00001) {SESSOUTV}
469
IF WLCMTD(1:1) .EQ. "A" THEN ! Method ANGLE
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)
477
REGRES/POLY {LINTAB} :AUX :XROT {DEGSTART}
478
PIXEL(1) = M$ABS(OUTPUTD(2)*wlcstep)
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)
485
REGRES/POLY {LINTAB} :AUX :X,:ORDER {DEGSTART},{DEGSTART}
486
PIXEL(1) = M$ABS(OUTPUTD(2)*wlcstep)
489
SAVE/REGRES {LINTAB} REGR
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))
496
IF {LWCL(1)} .GT. 0. THEN
497
PIXEL(2) = LWCL(1)* PIXEL(1)
499
PIXEL(2) = M$ABS(LWCL(1))
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
512
alpha(3) = outputr(1)
514
INQUIRE/KEY ANSW/C/1/3 "Start anyway (y/n, default no) ?"
515
IF ANSW(1:1) .EQ. "Y" THEN
518
WRITE/DESCR middummc.tbl TBLCONTR/I/4/1 0
519
WRITE/DESCR middummr.tbl TBLCONTR/I/4/1 0
527
! ************************ Start with a guess *************************
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:
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
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
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}:)}
557
WRITE/KEY GUETAB {GUESS}LINE.tbl
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)
567
PIXEL(2) = M$ABS(LWCL(1))
570
! Determines absolute order numbers in the line.tbl
571
! and updates the ORDPOS keyword.
573
define/local doris/I/1/1 0
574
doris = M$EXISTD(GUETAB,"RORDI")
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
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)
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
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"
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.
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}
626
WRITE/OUT "Cross-correlation shift value : 'SHIFT'"
628
! Applies the shift to the guess table and computes new coefficients
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
639
OUTPUTR(5) = {{LINTAB},REGRR(5)}
640
NBR(1) = {{LINTAB},REGRI(1)}
641
RMS(1) = 2.*OUTPUTR(5)/(ORDPOS(1)+ORDPOS(2))
643
WRITE/DESCR {LINTAB} ORDER/I/1/2 {ORDPOS(1)},{ORDPOS(2)}
644
WRITE/DESCR {LINTAB} PIXEL/R/1/2 {PIXEL(1)},{PIXEL(2)}
647
! ************************ Load table in overlay plane *****************
651
DEFINE/LOCAL NBROW/I/1/1 0 ? +lower
653
IF WLCVISU(1:1) .NE. "Y" RETURN
655
SELECT/TABLE &r (sequence.ne.0) {SESSOUTV}
657
COMPUTE/TABLE &r :SEQU = sequence
658
WRITE/DESCR middummr.tbl TBLCONTR/I/4/1 {NBROW}
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}
667
! ************************ Cursor interaction **************************
671
DEFINE/PARAMETER P1 HAY C "Method"
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)}"
679
IF WLCVISU(1:1) .NE. "Y" THEN
686
! If the command CENTER/MOMENT is modified here, check consistency with
687
! command used near label RECLICK
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"
695
CENTER/MOMENT CURSOR &c
697
COMPUTE/TABLE middummc :XERR = ABS(:XEND - :XSTART)/2.
698
COMPUTE/TABLE middummc :YERR = ABS(:YEND - :YSTART)/2.
700
SELECT/TABLE middummc (ABS(:STATUS).LT.0.01) {SESSOUTV}
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}
707
WRITE/DESCR middummr.tbl WLCMTD/C/1/8 {METHOD(1:5)}
711
! ********************* Prepares association loop values *******************
715
COMPUTE/TABLE {LINTAB} :IDENT = 0.
717
LAST(1) = {middummr.tbl,TBLCONTR(4)}
721
! ******************** Associates positions ****************************
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}}
732
select/table {LINTAB} ABS(:YNEW-'xypos(2)').LT.'xymax(2)'.and.ABS(:X-'xypos(1)').LT.'xymax(1)' {SESSOUTV}
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}
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
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)
765
! ******************************************************************
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.
777
WRITE/OUT "Iterative Identification:"
781
COMPUTE/REGRES {LINTAB} :WAVEC = REGR
782
IF WLCMTD(1:1) .NE. "C" THEN
783
COMPUTE/TABLE {LINTAB} :WAVEC = :WAVEC/:ORDER
785
COMPUTE/TABLE {LINTAB} :IDENT = 0.
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
796
WRITE/KEY PROGSTAT/I/1/1 0
798
IF PROGSTAT(1) .NE. 0 GOTO EXIT ! Problems in matching loop
801
!KB SELECT/TAB {LINTAB} (:IDENT.NE.0.0) {SESSOUTV}
802
SELECT/TAB {LINTAB} (ABS(:IDENT).GT.0.00001) {SESSOUTV}
805
WRITE/OUT "*** Number of identifications: {OUTPUTI(1)}"
807
IF WLCVISU(1:1) .EQ. "Y" THEN
809
LOAD/TAB {LINTAB} :X :YNEW {SESSOUTV}
812
STAT/TABLE {LINTAB} :RESIDUAL {SESSOUTV}
814
PIXEL(3) = M$ABS(RMS(2)/PIXEL(1))
816
! Tests whether one stops the iteration
817
IF ITER .EQ. NITER GOTO ENDIT
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
828
WRITE/OUT "*** STD. DEV. {RMS(2)} wavelength units, i.e. {PIXEL(3)} pixels ***"
830
IF PIXEL(3) .GT. LWCL(3) THEN
831
WRITE/OUT "Sorry, wrong identifications..."
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}
841
IF NBR(2) .GT. NBRMIN THEN
842
WRITE/KEY WLCMTD TWO-D
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}
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}
852
IF WLCMTD(1:1) .EQ. "C" REGRES/POLY {LINTAB} :IDENT :YNEW {DEGREE} KEYLONG {SESSOUTV}
854
SAVE/REGRES {LINTAB} REGR KEYLONG
863
!*****************************************************************************
867
SELECT/TABLE {LINTAB} (:SELECT.GT.1) {SESSOUTV}
868
WRITE/KEY IN_A {LINTAB}
869
WRITE/KEY IN_B {WLCMTD}
871
RUN STD_EXE:necwlcopt
872
IF OUTPUTI(1) .NE. 0 RETURN/EXIT
876
!************************************************************************
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."
884
WRITE/OUT "Please identify line {line} again."
887
CENTER/MOMENT CURSOR &c
889
COMPUTE/TABLE middummc :XERR = ABS(:XEND - :XSTART)/2.
890
COMPUTE/TABLE middummc :YERR = ABS(:YEND - :YSTART)/2.
892
SELECT/TABLE middummc (ABS(:STATUS).LT.0.01) {SESSOUTV}
896
DEFINE/LOCAL LINPAIR/I/1/1 0 ? +lower
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}
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}
923
!***********************************************************************
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.
933
convert = relord + absord
936
echord(2) = convert - echord(1)
937
echord(3) = convert - 1
939
ordpos(1) = echord(3)
940
ordpos(2) = echord(2)
942
copy/kd ordpos/I/1/2 {lintab} ORDER/I/1/2
944
compute/table {LINTAB} :ORDER = {CONVERT} - :Y
946
SELECT/TABLE {LINTAB} :Y.EQ.{relord} {SESSOUTV}
948
stat/table {LINTAB} :YNEW {SESSOUTV}
951
dispersion/hough ? ? {IMSIZE(1)},{{WLC},START(1)},{{WLC},STEP(1)} {LINTAB} {SESSOUTV}
953
COMPUTE/TABLE {LINTAB} :AUX = :WAVEC*:ORDER
954
SELECT/TABLE {LINTAB} :Y.EQ.{relord} {SESSOUTV}
955
REGRESS/POLY {LINTAB} :AUX :X 1 {SESSOUTV}
961
! First estimate on the error on the angle = 1. degree
962
factor = slop*M$TAN(1)
964
SAVE/REGR {LINTAB} REGR
966
COMPUTE/TABLE {LINTAB} :ERROR = (:YNEW-{YMEAN})*{FACTOR}
967
COMPUTE/TABLE {LINTAB} :ERROR = SQRT(:ERROR*:ERROR + {RMS})
968
COMPUTE/TABLE {LINTAB} :ERROR = :ERROR/:ORDER
970
stat/tab {lintab} :ERROR {SESSOUTV}
971
write/out "Error range: {outputr(1)} to {outputr(2)} wav. units"
973
SELECT/TABLE {LINTAB} ALL
978
! ******************************************************************
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.
988
WRITE/OUT ".. Iterative Identification:"
990
IF WLCVISU(1:1) .EQ. "Y" THEN
998
COMPUTE/REGRES {LINTAB} :WAVEC = REGR
999
COMPUTE/TABLE {LINTAB} :WAVEC = :WAVEC/:ORDER
1003
COMPUTE/TABLE {LINTAB} :IDENT = 0.
1004
COMPUTE/TABLE {LINTAB} :SELECT = 0.
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
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
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 ***"
1025
IF WLCVISU(1:1) .EQ. "Y" THEN
1027
PLOT/HISTO {LINTAB} :RESIDUAL
1028
select/table {LINTAB} :IDENT.GT.0.0
1029
PLOT/TABLE {LINTAB} :IDENT :RESIDUAL
1032
!KB SELECT/TAB {LINTAB} (:IDENT.NE.0.0) {SESSOUTV}
1033
SELECT/TAB {LINTAB} (ABS(:IDENT).GT.0.00001) {SESSOUTV}
1036
WRITE/OUT "*** Iteration: {ITER}. Number of identifications: {OUTPUTI(1)}"
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}
1043
LOAD/TAB {LINTAB} :X :YNEW {SESSOUTV}
1046
COMPUTE/TABLE {LINTAB} :AUX = :ORDER * :IDENT
1047
SELECT/TAB {LINTAB} (:IDENT.GT.0.0) {SESSOUTV}
1049
IF NBR(2) .GT. NBRMIN THEN
1050
WRITE/KEY WLCMTD TWO-D
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)
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}
1071
!Computes :WAVEC for all cases
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)}
1077
! Determines median of residuals
1079
select/table {lintab} :IDENT.GT.0.0
1081
IF WLCVISU(1:1) .EQ. "Y" THEN
1082
set/graph colour=1 pmode=1
1083
PLOT/TABLE {LINTAB} :IDENT :RESIDUAL
1086
copy/table {LINTAB} &l
1087
compute/table &l :R2 = ABS(:RESIDUAL)
1089
inputi(1) = {middumml.tbl,TBLCONTR(4)}/2
1090
outputr(4) = {middumml.tbl,:R2,@{inputi(1)}}
1091
RMS(2) = OUTPUTR(4)*2.5
1093
SELECT/TABLE {LINTAB} :IDENT.GT.0.0.AND.ABS(:RESIDUAL).LT.{RMS(2)} {SESSOUTV}
1095
IF WLCVISU(1:1) .EQ. "Y" THEN
1096
select/table {lintab} .NOT.SELECT
1098
OVER/TABLE {lintab} :IDENT :RESIDUAL
1099
select/table {lintab} .NOT.SELECT
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 ***"
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}
1119
REGRES/POLY {LINTAB} :AUX :X,:ORDER {DEGREE},{DEGREE} KEYLONG {SESSOUTV}
1120
SAVE/REGR {LINTAB} REGR KEYLONG {SESSOUTV}
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)}"
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
1139
IF PIXEL(3) .GT. LWCL(3) THEN
1140
WRITE/OUT "Sorry, wrong identifications..."
1153
!***************************************************************************
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
1162
SELECT/TABLE {LINTAB} ALL
1164
STATIST/IMA &h [<,@3:>,@3] OPTION=HN
1166
SELECT/TABLE {LINTAB} :PEAK.GT.{OUTPUTR(8)} {SESSOUTV} ! OUTPUTR(8) = Median
1169
SELECT/TABLE {LINTAB} SELECT.AND.:SELECT.GE.1 {SESSOUTV}
1172
SELECT/TABLE {LINTAB} ALL
1174
RATIO = 100.*SELECT/BRIGHT
1178
WRITE/OUT "************************** Verification ***********************************"
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."
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."
1193
!**********************************************************************
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
1201
select/table {LINTAB} all
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
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
1213
WRITE/KEY INPUTI/I/1/2 1,-10
1217
compute/table {LINTAB} :ERROR = 2.*SQRT(:BKG)
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