← Back to branch summary

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

« back to all changes in this revision

Viewing changes to stdred/echelle/src/necexor.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/necexor.for
 
1
C @(#)necexor.for       19.1 (ESO-DMD) 02/25/03 14:20:22
 
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 Massachusetts Ave, Cambridge, 
 
18
C MA 02139, USA.
 
19
C
 
20
C Correspondence 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: Copyright (c) 1987 European Southern Observatory,
 
31
C                                         all rights reserved
 
32
C
 
33
C.VERSION: 1.0  ESO-FORTRAN Conversion, AA  21:13 - 3 DEC 1987
 
34
C
 
35
C.LANGUAGE: F77+ESOext
 
36
C
 
37
C.AUTHOR: D.PONZ
 
38
C
 
39
C.IDENTIFICATION
 
40
C  program  ECHEXOR  version 1.2 840114
 
41
C
 
42
C.KEYWORDS
 
43
C
 
44
C  ECHELLE, CASPEC, DATA EXTRACTION
 
45
C
 
46
C.PURPOSE
 
47
C
 
48
C  EXTRACT ECHELLE ORDERS IN CASPEC IMAGES. ORDER POSITION ARE DEFINED
 
49
C  BY COEFFICIENTS OF A REGRESSION. SLIT WIDTH, ANGLE AND OFFSET
 
50
C  ARE DEFINED BY PARAMETERS
 
51
C            WIDTH,ANGLE,OFFSET
 
52
C
 
53
C.ALGORITHM
 
54
C
 
55
C   ORDERS ARE EXTRACTED BY PASSING A NUMERICAL SLIT OF GIVEN WIDTH 
 
56
C  AND ANGLE IN POSITIONS DEFINED BY THE REGRESSION COEFFS. THE SLIT
 
57
C  WILL SAMPLE AT INCREMENTS OF ONE STEP IN X.
 
58
C
 
59
C  THERE ARE THREE METHODS :
 
60
C   LINEAR   - PIXEL VALUES ARE ADDED, LINEAR INTERPOLATION
 
61
C   AVERAGE  - PIXEL VALUES ARE AVERAGED, LINEAR INTERPOLATION
 
62
C   WEIGHTED - PIXEL VALUES ARE WEIGHTD PROPORTIONALLY TO THE
 
63
C              PROFI PERPENDICULAR TO THE DISPERSION DIRECTION
 
64
C.INPUT/OUTPUT
 
65
C
 
66
C   EXTRACT/ORDER output = input width,angle,offset 
 
67
C                                table coeffs [ord1,ord2]
 
68
C  the connection file contains
 
69
C  IN_A   input FRM
 
70
C  IN_B   input table with coeffs
 
71
C  OUT_A  output FRM
 
72
C  INPUTR  width,angle
 
73
C  INPUTC  coeffs
 
74
C  INPUTI  ord1,ord2
 
75
C  ACTION               T - table option
 
76
C                       B - batch order by order
 
77
C 
 
78
C 
 
79
C 010123                last modif
 
80
C 
 
81
C------------------------------------------------------------------
 
82
C
 
83
      PROGRAM ECHEXO
 
84
C 
 
85
      IMPLICIT NONE
 
86
C 
 
87
      INTEGER NDIM
 
88
      PARAMETER (NDIM=2)
 
89
      INTEGER MADRID, NAC, NAR, IMNO
 
90
      INTEGER NWDTH
 
91
      INTEGER I,IAV,INDX,IORD,IORD1
 
92
      INTEGER IORD2,ISTAT,NA,NAXISA,NAXISB
 
93
      INTEGER NCOL,NPAR,NROW,NSORT
 
94
      INTEGER*8 PNTRA
 
95
      INTEGER*8 PNTRB
 
96
      INTEGER*8 PNTRC
 
97
      INTEGER   STATUS,TID,INO,TMPID
 
98
      INTEGER NPIXA(3),NPIXB(3),IPAR(10),KUN,KNUL
 
99
C
 
100
      DOUBLE PRECISION STEPA(3),STEPB(3), STARTA(3),STARTB(3)
 
101
      REAL  RPAR(10)
 
102
      REAL  CUTS(6)
 
103
      REAL XMIN,XMAX,XMIN1,XMAX1
 
104
      REAL PARAMS(12)
 
105
C
 
106
      DOUBLE PRECISION DPAR(100)
 
107
C
 
108
      CHARACTER*72   IDENA,IDENB
 
109
      CHARACTER*64   CUNITA,CUNITB,DESNAM
 
110
      CHARACTER*1    METH
 
111
      CHARACTER*60   FRMIN,FRMOUT,INA1,INA2
 
112
      CHARACTER*60   TABLE,NAME,TYPE
 
113
C
 
114
      INCLUDE 'MID_INCLUDE:ST_DEF.INC'
 
115
      COMMON /VMR/MADRID(1)
 
116
      INCLUDE 'MID_INCLUDE:ST_DAT.INC'
 
117
C
 
118
      DATA CUNITB/
 
119
     +     'FLUX            PIXEL                           ORDER'/
 
120
C
 
121
C  initialize system
 
122
C
 
123
      CALL STSPRO('ECHEXO')
 
124
      CALL STKRDC('IN_A',1,1,60,I,INA1,KUN,KNUL,STATUS)
 
125
      CALL STKRDC('IN_B',1,1,60,I,TABLE,KUN,KNUL,STATUS)
 
126
      CALL STKRDC('IN_C',1,1,1,I,METH,KUN,KNUL,STATUS)
 
127
      CALL STKRDC('OUT_A',1,1,60,I,INA2,KUN,KNUL,STATUS)
 
128
      CALL STKRDC('INPUTC',1,1,60,I,NAME,KUN,KNUL,STATUS)
 
129
      CALL STKRDR('INPUTR',1,3,I,RPAR,KUN,KNUL,STATUS)
 
130
      CALL STKRDI('INPUTI',1,2,I,IPAR,KUN,KNUL,STATUS)
 
131
      IORD1  = IPAR(1)
 
132
      IORD2  = IPAR(2)
 
133
      FRMIN = INA1
 
134
      FRMOUT = INA2
 
135
      IF (METH(1:1) .EQ. 'l') METH = 'L'
 
136
      IF (METH(1:1) .EQ. 'a') METH = 'A'
 
137
      IF (METH(1:1) .EQ. 'm') METH = 'M'
 
138
C
 
139
C  input parameters in table
 
140
C
 
141
      CALL TBTOPN(TABLE,F_U_MODE,TID,STATUS)
 
142
      CALL TBIGET(TID,NCOL,NROW,NSORT,NAC,NAR,STATUS)
 
143
C
 
144
C  map input FRM
 
145
C
 
146
      CALL STIGET(FRMIN,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE,
 
147
     .            NDIM,NAXISA,NPIXA,STARTA,STEPA,IDENA,
 
148
     +            CUNITA,PNTRA,IMNO,STATUS)
 
149
      PARAMS(1) = RPAR(1)/ABS(STEPA(2))
 
150
      PARAMS(2) = RPAR(2)
 
151
      PARAMS(3) = RPAR(3)/ABS(STEPA(2))
 
152
 
 
153
      IF (METH(1:1) .EQ. 'M') THEN
 
154
         NWDTH = PARAMS(1) + 2
 
155
         CALL STFCRE('column',D_R4_FORMAT,F_X_MODE,1,NWDTH,TMPID,STATUS)
 
156
         CALL STFMAP(TMPID,F_X_MODE,1,NWDTH,IAV,PNTRC,STATUS)
 
157
      ENDIF
 
158
C
 
159
C  read coeffs
 
160
C
 
161
      NAME   = NAME(1:59)//' '
 
162
      INDX   = INDEX(NAME,' ') - 1
 
163
      DESNAM = NAME(1:INDX)//'C'
 
164
      CALL STDRDC(TID,DESNAM,1,17,4,IAV,TYPE,KUN,KNUL,ISTAT)
 
165
      TYPE   = 'MULT'
 
166
      IF (TYPE.EQ.'MULT') THEN
 
167
         DESNAM = NAME(1:INDX)//'I'
 
168
         CALL STDRDI(TID,DESNAM,4,4,IAV,IPAR,KUN,KNUL,ISTAT)
 
169
         NA     = (IPAR(3)+1)* (IPAR(4)+1)
 
170
         DESNAM = NAME(1:INDX)//'R'
 
171
         CALL STDRDR(TID,DESNAM,1,4,IAV,PARAMS(4),KUN,KNUL,ISTAT)
 
172
         DESNAM = NAME(1:INDX)//'D'
 
173
         CALL STDRDD(TID,DESNAM,1,NA,IAV,DPAR,KUN,KNUL,ISTAT)
 
174
         CALL TBTCLO(TID,STATUS)
 
175
C
 
176
C ... map output FRM
 
177
C
 
178
         IF (IORD1.EQ.IORD2 .AND. IORD1.EQ.0) THEN
 
179
            IORD1  = PARAMS(6)
 
180
            IORD2  = PARAMS(7)
 
181
         END IF
 
182
         STARTB(1) = STARTA(1)
 
183
         STARTB(2) = IORD1
 
184
         STEPB(1)  = STEPA(1)
 
185
         STEPB(2)  = 1.
 
186
         IDENB     = IDENA
 
187
         NAXISB    = 2
 
188
         NPIXB(1)  = NPIXA(1)
 
189
         NPIXB(2)  = IORD2 - IORD1 + 1
 
190
         CALL STIPUT(FRMOUT,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE,
 
191
     +        NAXISB,NPIXB,STARTB,
 
192
     +        STEPB,IDENB,CUNITB,PNTRB,INO,STATUS) 
 
193
C     
 
194
C  extract order
 
195
C
 
196
         NPAR      = 10
 
197
         PARAMS(8) = IPAR(3)
 
198
         PARAMS(9) = IPAR(4)
 
199
         XMIN      = 1.E30
 
200
         XMAX      = -XMIN
 
201
         DO 10 I = IORD1,IORD2
 
202
            
 
203
            PARAMS(10) = I
 
204
            IORD       = I - IORD1 + 1
 
205
            
 
206
            IF (METH(1:1) .EQ. 'M') THEN
 
207
               CALL EXTR2M(MADRID(PNTRA),NPIXA(1),NPIXA(2),
 
208
     +              MADRID(PNTRB),NPIXB(1),NPIXB(2),MADRID(PNTRC),NWDTH,
 
209
     +              IORD,NPAR,PARAMS,NA,DPAR,XMIN1,XMAX1,METH,
 
210
     +              STARTA,STEPA)
 
211
            ELSE
 
212
               CALL EXTR1M(MADRID(PNTRA),NPIXA(1),NPIXA(2),
 
213
     +              MADRID(PNTRB),NPIXB(1),NPIXB(2),
 
214
     +              IORD,NPAR,PARAMS,NA,DPAR,XMIN1,XMAX1,METH,
 
215
     +              STARTA,STEPA)
 
216
            ENDIF
 
217
 
 
218
            XMIN   = AMIN1(XMIN,XMIN1)
 
219
            XMAX   = AMAX1(XMAX,XMAX1)
 
220
 10      CONTINUE
 
221
      ELSE
 
222
         CALL STTPUT(' Error in dispersion coefficients ')
 
223
         GO TO 20
 
224
      END IF
 
225
C
 
226
C  descriptor handling
 
227
C
 
228
      CUTS(1) = XMIN
 
229
      CUTS(2) = XMAX
 
230
      CUTS(3) = XMIN
 
231
      CUTS(4) = XMAX
 
232
      CUTS(5) = 0.
 
233
      CUTS(6) = 0.
 
234
      CALL DSCUPT(IMNO,INO,' ',STATUS)
 
235
      CALL STDWRR(INO,'LHCUTS',CUTS,1,6,KUN,STATUS)
 
236
C
 
237
C  free data
 
238
C
 
239
C      CALL STFCLO(INO,STATUS)
 
240
 20   CALL STSEPI
 
241
      END
 
242
 
 
243
      SUBROUTINE EXTR1M(IN,NPIXA1,NPIXA2,OUT,NB1,NB2,INDX,NPAR,PAR,ND,
 
244
     +     DPAR,XMIN,XMAX,METH,STARTA,STEPA)
 
245
C
 
246
C       EXTRACTION ROUTINE - BILINEAR INTERPOLATION
 
247
C
 
248
C   ARGUMENTS
 
249
C IN     REAL*4  INPUT IMAGE
 
250
C NPIXA1 INTG*4  NO. OF SAMPLES
 
251
C NPIXA2 INTG*4  NO. OF LINES
 
252
C OUT    REAL*4  EXTRACTED ORDERS
 
253
C NB1    INTG*4  NO. OF SAMPLES
 
254
C NB2    INTG*4  NO. OF ORDERS
 
255
C INDX   INTG*4  SEQUENTIAL NO. OF THE ORDER
 
256
C NPAR   INTG*4  NO. OF EXTRACTION PARAMETERS
 
257
C PAR    REAL*4  EXTRACTION PARAMETERS
 
258
C ND     INTG*4  NO. OF COEFFS.
 
259
C DPAR   REAL*8  COEFFS. TO DEFINE ORDER POSITION
 
260
C XMIN   REAL*4  MINIMUM VALUE
 
261
C XMAX   REAL*4  MAXIMUMVALUE
 
262
C
 
263
      IMPLICIT NONE
 
264
C 
 
265
      INTEGER  NB1, NB2, NPAR, NPIXA1, NPIXA2, ND
 
266
      INTEGER  LOSLIT, HISLIT
 
267
      INTEGER  NWW, IWIDTH, IYLIM, NX, N1, N2, K, L, IX, INDX
 
268
      INTEGER  IX1, IP, L1, K1, IW
 
269
C 
 
270
      REAL OUT(NB1,NB2),PAR(NPAR)
 
271
      REAL IN(NPIXA1,NPIXA2)
 
272
      REAL XMIN,XMAX
 
273
C 
 
274
      DOUBLE PRECISION DPAR(ND),STARTA(2),STEPA(2)
 
275
      DOUBLE PRECISION XVAL,YVAL,AX
 
276
      DOUBLE PRECISION WIDTH,SKEW,HWIDTH,OFFSET,DC
 
277
      DOUBLE PRECISION YV,WW(50)
 
278
      DOUBLE PRECISION SKWDX,SKWDY,SKWOF,SUM
 
279
      DOUBLE PRECISION DEGREE,DD,DSIN,DCOS
 
280
      DOUBLE PRECISION YLOW, YHIGH, WSUM
 
281
C 
 
282
      LOGICAL IF1D
 
283
C 
 
284
      CHARACTER*1      METH
 
285
C 
 
286
      DOUBLE PRECISION IX_R8R8, XI_R8R8
 
287
C
 
288
C     DEGREE = RAD / DEGREE
 
289
C
 
290
      DATA DEGREE/0.017453292519943295769237D+0/
 
291
C
 
292
 
 
293
      NWW    = 50
 
294
C
 
295
C ... define extraction constants
 
296
C
 
297
      WIDTH  = PAR(1)
 
298
      IWIDTH = WIDTH
 
299
      HWIDTH = 0.5*WIDTH
 
300
      SKEW   = PAR(2)
 
301
      SKWDX  = 0.0
 
302
      SKWDY  = 1.0
 
303
      SKWOF  = 0.0
 
304
      XMIN   = 0.0
 
305
      XMAX   = XMIN
 
306
C      OFFSET = PAR(3) - HWIDTH + 0.5
 
307
      OFFSET = PAR(3)
 
308
      IF (SKEW.NE.0.0) THEN       
 
309
         DD      = SKEW*DEGREE
 
310
         SKWDX  = -DSIN(DD)     !   DX PER PIXEL
 
311
         SKWDY  = DCOS(DD)      !   DY PER PIXEL
 
312
         SKWOF  = -SKWDX*HWIDTH !   2D INTERPOLATION
 
313
         IF1D   = .FALSE.
 
314
      END IF
 
315
      AX     = 1.D0 + SKWOF
 
316
      IYLIM  = NPIXA2 - IWIDTH
 
317
      NX     = SKWOF
 
318
      N1     = NX
 
319
      N2     = NPIXA1 - NX
 
320
C
 
321
C ... get ordinate values
 
322
C
 
323
      K      = PAR(8)
 
324
      L      = PAR(9)
 
325
      YV     = PAR(10)
 
326
C
 
327
C ... loop over the samples
 
328
C
 
329
      DO 10 IX = 1,N1
 
330
         OUT(IX,INDX) = 0.
 
331
 10   CONTINUE
 
332
 
 
333
      DO 100 IX = N2 + 1,NPIXA1
 
334
         OUT(IX,INDX) = 0.
 
335
 100  CONTINUE
 
336
 
 
337
      DO 50 IX1 = N1 + 1,N2
 
338
 
 
339
         IP     = 0
 
340
         DC     = 1.D0
 
341
         YVAL   = 0.D0
 
342
         XVAL   = IX_R8R8(AX,STARTA(1),STEPA(1))
 
343
 
 
344
C Computes position of the center of the slit from the coefficients of
 
345
C the bivariate polynomial stored in table order.
 
346
 
 
347
         DO 30 L1 = 0,L
 
348
            IP     = IP + 1
 
349
            WW(IP) = DC
 
350
            YVAL   = YVAL + WW(IP)*DPAR(IP)
 
351
            DO 20 K1 = 1,K
 
352
               IP     = IP + 1
 
353
               WW(IP) = WW(IP-1)*XVAL
 
354
               YVAL   = YVAL + WW(IP)*DPAR(IP)
 
355
 20         CONTINUE
 
356
            DC     = DC*YV
 
357
 30      CONTINUE
 
358
C
 
359
C Corrects the central position for the offset
 
360
C OFFSET does not take into account HWIDTH 
 
361
C
 
362
         YVAL   = XI_R8R8(YVAL,STARTA(2),STEPA(2)) + OFFSET
 
363
C     
 
364
         YLOW   = YVAL - HWIDTH
 
365
         YHIGH  = YVAL + HWIDTH
 
366
         LOSLIT = INT(YLOW)
 
367
         HISLIT = INT(YHIGH)
 
368
 
 
369
         IF (LOSLIT.LT.1)      LOSLIT = 1
 
370
         IF (HISLIT.GT.NPIXA2) HISLIT = NPIXA2
 
371
 
 
372
         SUM    = 0.D0
 
373
         WSUM   = 0.D0
 
374
         IX     = INT(AX)
 
375
C
 
376
         IF (LOSLIT.EQ.HISLIT) THEN
 
377
            SUM = (YHIGH-YLOW)*IN(IX,LOSLIT)
 
378
C
 
379
         ELSE
 
380
C
 
381
            DO 40 IW = LOSLIT,HISLIT
 
382
               
 
383
               IF (IW.EQ.LOSLIT) THEN
 
384
                  SUM = SUM + IN(IX,IW)*(LOSLIT+1.-YLOW)
 
385
                  
 
386
               ELSE IF (IW.EQ.HISLIT) THEN
 
387
                  SUM = SUM + IN(IX,IW)*(YHIGH-HISLIT)
 
388
 
 
389
               ELSE
 
390
                  SUM = SUM + IN(IX,IW)
 
391
 
 
392
               ENDIF
 
393
 
 
394
 
 
395
 40         CONTINUE
 
396
C
 
397
         ENDIF
 
398
C
 
399
         IF (METH(1:1) .EQ. 'L') THEN ! Linear extraction
 
400
            OUT(IX1,INDX) = SUM
 
401
         ELSE                   ! Average extraction
 
402
            OUT(IX1,INDX) = SUM/WIDTH
 
403
         ENDIF
 
404
         XMIN   = AMIN1(XMIN,OUT(IX1,INDX))
 
405
         XMAX   = AMAX1(XMAX,OUT(IX1,INDX))
 
406
         AX     = AX + 1.D0
 
407
 50   CONTINUE
 
408
 
 
409
      RETURN
 
410
 
 
411
      END
 
412
 
 
413
 
 
414
 
 
415
      SUBROUTINE EXTR2M(IN,NPIXA1,NPIXA2,OUT,NB1,NB2,COL,NWDTH,INDX,NPAR
 
416
     $     ,PAR,ND,DPAR,XMIN,XMAX,METH,STARTA,STEPA)
 
417
C
 
418
C       EXTRACTION ROUTINE - MEDIAN
 
419
C
 
420
C   ARGUMENTS
 
421
C IN     REAL*4  INPUT IMAGE
 
422
C NPIXA1 INTG*4  NO. OF SAMPLES
 
423
C NPIXA2 INTG*4  NO. OF LINES
 
424
C OUT    REAL*4  EXTRACTED ORDERS
 
425
C NB1    INTG*4  NO. OF SAMPLES
 
426
C NB2    INTG*4  NO. OF ORDERS
 
427
C INDX   INTG*4  SEQUENTIAL NO. OF THE ORDER
 
428
C NPAR   INTG*4  NO. OF EXTRACTION PARAMETERS
 
429
C PAR    REAL*4  EXTRACTION PARAMETERS
 
430
C ND     INTG*4  NO. OF COEFFS.
 
431
C DPAR   REAL*8  COEFFS. TO DEFINE ORDER POSITION
 
432
C XMIN   REAL*4  MINIMUM VALUE
 
433
C XMAX   REAL*4  MAXIMUMVALUE
 
434
C
 
435
      IMPLICIT NONE
 
436
C 
 
437
      INTEGER  NB1, NB2, NPAR, NPIXA1, NPIXA2, NWDTH, NWDTH2, ND
 
438
      INTEGER  NWW, IWIDTH, IYLIM, NX, N1, N2, K, L, IX, INDX
 
439
      INTEGER  IX1, IP, L1, K1, IW
 
440
      INTEGER  LOSLIT, HISLIT
 
441
C 
 
442
      REAL OUT(NB1,NB2),PAR(NPAR)
 
443
      REAL IN(NPIXA1,NPIXA2),COL(NWDTH)
 
444
      REAL XMIN,XMAX,LIMITS(2)
 
445
C 
 
446
      DOUBLE PRECISION DPAR(ND),STARTA(2),STEPA(2)
 
447
      DOUBLE PRECISION XVAL,YVAL,AX
 
448
      DOUBLE PRECISION WIDTH,SKEW,HWIDTH,OFFSET,DC
 
449
      DOUBLE PRECISION YV,WW(50)
 
450
      DOUBLE PRECISION SKWDX,SKWDY,SKWOF,SUM
 
451
      DOUBLE PRECISION DEGREE,DD,DSIN,DCOS
 
452
      DOUBLE PRECISION YLOW, YHIGH, WSUM
 
453
      DOUBLE PRECISION IX_R8R8, XI_R8R8
 
454
C 
 
455
      LOGICAL IF1D
 
456
C 
 
457
      CHARACTER*1      METH
 
458
      CHARACTER*80   OUTPUT
 
459
C
 
460
C     DEGREE = RAD / DEGREE
 
461
C
 
462
      DATA DEGREE/0.017453292519943295769237D+0/
 
463
      DATA LIMITS/0.0,0.0/
 
464
C     
 
465
 
 
466
      NWW    = 50
 
467
C
 
468
C ... define extraction constants
 
469
C
 
470
      WIDTH  = PAR(1)
 
471
      IWIDTH = WIDTH
 
472
      HWIDTH = 0.5*WIDTH
 
473
      NWDTH2 = HWIDTH
 
474
      SKEW   = PAR(2)
 
475
      SKWDX  = 0.0
 
476
      SKWDY  = 1.0
 
477
      SKWOF  = 0.0
 
478
      XMIN   = 0.0
 
479
      XMAX   = XMIN
 
480
C      OFFSET = PAR(3) - HWIDTH + 0.5
 
481
      OFFSET = PAR(3)
 
482
      IF (SKEW.NE.0.0) THEN       
 
483
          DD      = SKEW*DEGREE
 
484
          SKWDX  = -DSIN(DD)      !   DX PER PIXEL
 
485
          SKWDY  = DCOS(DD)       !   DY PER PIXEL
 
486
          SKWOF  = -SKWDX*HWIDTH  !   2D INTERPOLATION
 
487
          IF1D   = .FALSE.
 
488
      END IF
 
489
      AX     = 1.D0 + SKWOF
 
490
      IYLIM  = NPIXA2 - IWIDTH
 
491
      NX     = SKWOF
 
492
      N1     = NX
 
493
      N2     = NPIXA1 - NX
 
494
C
 
495
C ... get ordinate values
 
496
C
 
497
      K      = PAR(8)
 
498
      L      = PAR(9)
 
499
      YV     = PAR(10)
 
500
C
 
501
C ... loop over the samples
 
502
C
 
503
      DO 10 IX = 1,N1
 
504
         OUT(IX,INDX) = 0.
 
505
 10   CONTINUE
 
506
 
 
507
      DO 100 IX = N2 + 1,NPIXA1
 
508
         OUT(IX,INDX) = 0.
 
509
 100  CONTINUE
 
510
 
 
511
      DO 50 IX1 = N1 + 1,N2
 
512
 
 
513
         IP     = 0
 
514
         DC     = 1.D0
 
515
         YVAL   = 0.D0
 
516
         XVAL   = IX_R8R8(AX,STARTA(1),STEPA(1))
 
517
 
 
518
C Computes position of the center of the slit from the coefficients of
 
519
C the bivariate polynomial stored in table order.
 
520
 
 
521
         DO 30 L1 = 0,L
 
522
            IP     = IP + 1
 
523
            WW(IP) = DC
 
524
            YVAL   = YVAL + WW(IP)*DPAR(IP)
 
525
            DO 20 K1 = 1,K
 
526
               IP     = IP + 1
 
527
               WW(IP) = WW(IP-1)*XVAL
 
528
               YVAL   = YVAL + WW(IP)*DPAR(IP)
 
529
 20         CONTINUE
 
530
            DC     = DC*YV
 
531
 30      CONTINUE
 
532
C
 
533
C Corrects the central position for the offset
 
534
C OFFSET does not take into account HWIDTH 
 
535
C
 
536
         YVAL   = XI_R8R8(YVAL,STARTA(2),STEPA(2)) + OFFSET
 
537
C     
 
538
         YLOW   = YVAL - HWIDTH
 
539
         YHIGH  = YVAL + HWIDTH
 
540
         LOSLIT = INT(YLOW)
 
541
         HISLIT = INT(YHIGH)
 
542
 
 
543
         IF (LOSLIT.LT.1)      LOSLIT = 1
 
544
         IF (HISLIT.GT.NPIXA2) HISLIT = NPIXA2
 
545
 
 
546
         SUM    = 0.D0
 
547
         WSUM   = 0.D0
 
548
         IX     = INT(AX)
 
549
C
 
550
         IF (LOSLIT.EQ.HISLIT) THEN
 
551
            OUT(IX1,INDX) = IN(IX,LOSLIT)
 
552
         ELSE
 
553
            K1 = 0
 
554
            DO 40 IW = LOSLIT,HISLIT
 
555
               K1 = K1 + 1
 
556
               IF (IW.EQ.LOSLIT) THEN
 
557
                  COL(K1) = IN(IX,IW)!*(LOSLIT+1.-YLOW)
 
558
               ELSE IF (IW.EQ.HISLIT) THEN
 
559
                  COL(K1) = IN(IX,IW)!*MIN(1.0,YHIGH-HISLIT)
 
560
               ELSE
 
561
                  COL(K1) = IN(IX,IW)
 
562
               ENDIF
 
563
 40         CONTINUE            !determine MEDIAN:
 
564
            L1 = K1 / 2
 
565
            IF (K1.GT.0) THEN
 
566
               CALL GET_MEDIAN(COL,NWDTH,K1,L1,OUT(IX1,INDX),IW)
 
567
               IF (IW.LT.0) THEN
 
568
                  WRITE(OUTPUT,12345) IX1,INDX
 
569
                  CALL STTPUT(OUTPUT)
 
570
               ENDIF
 
571
            ELSE
 
572
               OUT(IX1,INDX) = -1.0
 
573
            ENDIF
 
574
         ENDIF
 
575
C
 
576
         XMIN   = AMIN1(XMIN,OUT(IX1,INDX))
 
577
         XMAX   = AMAX1(XMAX,OUT(IX1,INDX))
 
578
         AX     = AX + 1.D0
 
579
 50   CONTINUE
 
580
C
 
581
      RETURN
 
582
C 
 
583
12345 FORMAT('ERROR at ',I5,',',I5,' occurred')
 
584
C
 
585
      END
 
586
 
 
587
      SUBROUTINE GET_MEDIAN(RA,N,NUSE,LL,VALU,STATUS)
 
588
 
 
589
      IMPLICIT NONE
 
590
C 
 
591
      INTEGER   N,LL,STATUS
 
592
      INTEGER   NUSE
 
593
C 
 
594
      REAL      RA(N),VALU
 
595
C 
 
596
      VALU = 0.0
 
597
      STATUS = 0
 
598
C 
 
599
      IF (NUSE.LE.3) THEN
 
600
         IF (NUSE.LT.1) THEN
 
601
            STATUS = -1
 
602
            RETURN
 
603
         ELSE IF (NUSE.LE.2) THEN
 
604
            STATUS = 2
 
605
            VALU = RA(1)
 
606
c            VALU = (RA(1) + RA(2))/2.0               !for ndim = 1,2, median is 1st elem.
 
607
            RETURN
 
608
         ELSE 
 
609
            STATUS = 3
 
610
            VALU = RA(2)               !for ndim = 3, median is 2nd elem.
 
611
            RETURN
 
612
         ENDIF
 
613
      ENDIF
 
614
C 
 
615
      CALL NSORT(RA,NUSE,LL,VALU)
 
616
      RETURN
 
617
C 
 
618
      END