C @(#)necexor.for 17.1.1.1 (ESO-DMD) 01/25/02 17:51:29 C=========================================================================== C Copyright (C) 1995 European Southern Observatory (ESO) C C This program is free software; you can redistribute it and/or C modify it under the terms of the GNU General Public License as C published by the Free Software Foundation; either version 2 of C the License, or (at your option) any later version. C C This program is distributed in the hope that it will be useful, C but WITHOUT ANY WARRANTY; without even the implied warranty of C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the C GNU General Public License for more details. C C You should have received a copy of the GNU General Public C License along with this program; if not, write to the Free C Software Foundation, Inc., 675 Massachusetts Ave, Cambridge, C MA 02139, USA. C C Correspondence concerning ESO-MIDAS should be addressed as follows: C Internet e-mail: midas@eso.org C Postal address: European Southern Observatory C Data Management Division C Karl-Schwarzschild-Strasse 2 C D 85748 Garching bei Muenchen C GERMANY C=========================================================================== C C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C.COPYRIGHT: Copyright (c) 1987 European Southern Observatory, C all rights reserved C C.VERSION: 1.0 ESO-FORTRAN Conversion, AA 21:13 - 3 DEC 1987 C C.LANGUAGE: F77+ESOext C C.AUTHOR: D.PONZ C C.IDENTIFICATION C program ECHEXOR version 1.2 840114 C C.KEYWORDS C C ECHELLE, CASPEC, DATA EXTRACTION C C.PURPOSE C C EXTRACT ECHELLE ORDERS IN CASPEC IMAGES. ORDER POSITION ARE DEFINED C BY COEFFICIENTS OF A REGRESSION. SLIT WIDTH, ANGLE AND OFFSET C ARE DEFINED BY PARAMETERS C WIDTH,ANGLE,OFFSET C C.ALGORITHM C C ORDERS ARE EXTRACTED BY PASSING A NUMERICAL SLIT OF GIVEN WIDTH C AND ANGLE IN POSITIONS DEFINED BY THE REGRESSION COEFFS. THE SLIT C WILL SAMPLE AT INCREMENTS OF ONE STEP IN X. C C THERE ARE THREE METHODS : C LINEAR - PIXEL VALUES ARE ADDED, LINEAR INTERPOLATION C AVERAGE - PIXEL VALUES ARE AVERAGED, LINEAR INTERPOLATION C WEIGHTED - PIXEL VALUES ARE WEIGHTD PROPORTIONALLY TO THE C PROFI PERPENDICULAR TO THE DISPERSION DIRECTION C.INPUT/OUTPUT C C EXTRACT/ORDER output = input width,angle,offset C table coeffs [ord1,ord2] C the connection file contains C IN_A input FRM C IN_B input table with coeffs C OUT_A output FRM C INPUTR width,angle C INPUTC coeffs C INPUTI ord1,ord2 C ACTION T - table option C B - batch order by order C C C 010123 last modif C C------------------------------------------------------------------ C PROGRAM ECHEXO C IMPLICIT NONE C INTEGER NDIM PARAMETER (NDIM=2) INTEGER MADRID, NAC, NAR, IMNO INTEGER NWDTH INTEGER I,IAV,INDX,IORD,IORD1 INTEGER IORD2,ISTAT,NA,NAXISA,NAXISB INTEGER NCOL,NPAR,NROW,NSORT INTEGER*8 PNTRA INTEGER*8 PNTRB INTEGER*8 PNTRC INTEGER STATUS,TID,INO,TMPID INTEGER NPIXA(3),NPIXB(3),IPAR(10),KUN,KNUL C DOUBLE PRECISION STEPA(3),STEPB(3), STARTA(3),STARTB(3) REAL RPAR(10) REAL CUTS(6) REAL XMIN,XMAX,XMIN1,XMAX1 REAL PARAMS(12) C DOUBLE PRECISION DPAR(100) C CHARACTER*72 IDENA,IDENB CHARACTER*64 CUNITA,CUNITB,DESNAM CHARACTER*1 METH CHARACTER*60 FRMIN,FRMOUT,INA1,INA2 CHARACTER*60 TABLE,NAME,TYPE C INCLUDE 'MID_INCLUDE:ST_DEF.INC' COMMON /VMR/MADRID(1) INCLUDE 'MID_INCLUDE:ST_DAT.INC' C DATA CUNITB/ + 'FLUX PIXEL ORDER'/ C C initialize system C CALL STSPRO('ECHEXO') CALL STKRDC('IN_A',1,1,60,I,INA1,KUN,KNUL,STATUS) CALL STKRDC('IN_B',1,1,60,I,TABLE,KUN,KNUL,STATUS) CALL STKRDC('IN_C',1,1,1,I,METH,KUN,KNUL,STATUS) CALL STKRDC('OUT_A',1,1,60,I,INA2,KUN,KNUL,STATUS) CALL STKRDC('INPUTC',1,1,60,I,NAME,KUN,KNUL,STATUS) CALL STKRDR('INPUTR',1,3,I,RPAR,KUN,KNUL,STATUS) CALL STKRDI('INPUTI',1,2,I,IPAR,KUN,KNUL,STATUS) IORD1 = IPAR(1) IORD2 = IPAR(2) FRMIN = INA1 FRMOUT = INA2 IF (METH(1:1) .EQ. 'l') METH = 'L' IF (METH(1:1) .EQ. 'a') METH = 'A' IF (METH(1:1) .EQ. 'm') METH = 'M' C C input parameters in table C CALL TBTOPN(TABLE,F_U_MODE,TID,STATUS) CALL TBIGET(TID,NCOL,NROW,NSORT,NAC,NAR,STATUS) C C map input FRM C CALL STIGET(FRMIN,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE, . NDIM,NAXISA,NPIXA,STARTA,STEPA,IDENA, + CUNITA,PNTRA,IMNO,STATUS) PARAMS(1) = RPAR(1)/ABS(STEPA(2)) PARAMS(2) = RPAR(2) PARAMS(3) = RPAR(3)/ABS(STEPA(2)) IF (METH(1:1) .EQ. 'M') THEN NWDTH = PARAMS(1) + 2 CALL STFCRE('column',D_R4_FORMAT,F_X_MODE,1,NWDTH,TMPID,STATUS) CALL STFMAP(TMPID,F_X_MODE,1,NWDTH,IAV,PNTRC,STATUS) ENDIF C C read coeffs C NAME = NAME(1:59)//' ' INDX = INDEX(NAME,' ') - 1 DESNAM = NAME(1:INDX)//'C' CALL STDRDC(TID,DESNAM,1,17,4,IAV,TYPE,KUN,KNUL,ISTAT) TYPE = 'MULT' IF (TYPE.EQ.'MULT') THEN DESNAM = NAME(1:INDX)//'I' CALL STDRDI(TID,DESNAM,4,4,IAV,IPAR,KUN,KNUL,ISTAT) NA = (IPAR(3)+1)* (IPAR(4)+1) DESNAM = NAME(1:INDX)//'R' CALL STDRDR(TID,DESNAM,1,4,IAV,PARAMS(4),KUN,KNUL,ISTAT) DESNAM = NAME(1:INDX)//'D' CALL STDRDD(TID,DESNAM,1,NA,IAV,DPAR,KUN,KNUL,ISTAT) CALL TBTCLO(TID,STATUS) C C ... map output FRM C IF (IORD1.EQ.IORD2 .AND. IORD1.EQ.0) THEN IORD1 = PARAMS(6) IORD2 = PARAMS(7) END IF STARTB(1) = STARTA(1) STARTB(2) = IORD1 STEPB(1) = STEPA(1) STEPB(2) = 1. IDENB = IDENA NAXISB = 2 NPIXB(1) = NPIXA(1) NPIXB(2) = IORD2 - IORD1 + 1 CALL STIPUT(FRMOUT,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE, + NAXISB,NPIXB,STARTB, + STEPB,IDENB,CUNITB,PNTRB,INO,STATUS) C C extract order C NPAR = 10 PARAMS(8) = IPAR(3) PARAMS(9) = IPAR(4) XMIN = 1.E30 XMAX = -XMIN DO 10 I = IORD1,IORD2 PARAMS(10) = I IORD = I - IORD1 + 1 IF (METH(1:1) .EQ. 'M') THEN CALL EXTR2M(MADRID(PNTRA),NPIXA(1),NPIXA(2), + MADRID(PNTRB),NPIXB(1),NPIXB(2),MADRID(PNTRC),NWDTH, + IORD,NPAR,PARAMS,NA,DPAR,XMIN1,XMAX1,METH, + STARTA,STEPA) ELSE CALL EXTR1M(MADRID(PNTRA),NPIXA(1),NPIXA(2), + MADRID(PNTRB),NPIXB(1),NPIXB(2), + IORD,NPAR,PARAMS,NA,DPAR,XMIN1,XMAX1,METH, + STARTA,STEPA) ENDIF XMIN = AMIN1(XMIN,XMIN1) XMAX = AMAX1(XMAX,XMAX1) 10 CONTINUE ELSE CALL STTPUT(' Error in dispersion coefficients ') GO TO 20 END IF C C descriptor handling C CUTS(1) = XMIN CUTS(2) = XMAX CUTS(3) = XMIN CUTS(4) = XMAX CUTS(5) = 0. CUTS(6) = 0. CALL DSCUPT(IMNO,INO,' ',STATUS) CALL STDWRR(INO,'LHCUTS',CUTS,1,6,KUN,STATUS) C C free data C C CALL STFCLO(INO,STATUS) 20 CALL STSEPI END SUBROUTINE EXTR1M(IN,NPIXA1,NPIXA2,OUT,NB1,NB2,INDX,NPAR,PAR,ND, + DPAR,XMIN,XMAX,METH,STARTA,STEPA) C C EXTRACTION ROUTINE - BILINEAR INTERPOLATION C C ARGUMENTS C IN REAL*4 INPUT IMAGE C NPIXA1 INTG*4 NO. OF SAMPLES C NPIXA2 INTG*4 NO. OF LINES C OUT REAL*4 EXTRACTED ORDERS C NB1 INTG*4 NO. OF SAMPLES C NB2 INTG*4 NO. OF ORDERS C INDX INTG*4 SEQUENTIAL NO. OF THE ORDER C NPAR INTG*4 NO. OF EXTRACTION PARAMETERS C PAR REAL*4 EXTRACTION PARAMETERS C ND INTG*4 NO. OF COEFFS. C DPAR REAL*8 COEFFS. TO DEFINE ORDER POSITION C XMIN REAL*4 MINIMUM VALUE C XMAX REAL*4 MAXIMUMVALUE C IMPLICIT NONE C INTEGER NB1, NB2, NPAR, NPIXA1, NPIXA2, ND INTEGER LOSLIT, HISLIT INTEGER NWW, IWIDTH, IYLIM, NX, N1, N2, K, L, IX, INDX INTEGER IX1, IP, L1, K1, IW C REAL OUT(NB1,NB2),PAR(NPAR) REAL IN(NPIXA1,NPIXA2) REAL XMIN,XMAX C DOUBLE PRECISION DPAR(ND),STARTA(2),STEPA(2) DOUBLE PRECISION XVAL,YVAL,AX DOUBLE PRECISION WIDTH,SKEW,HWIDTH,OFFSET,DC DOUBLE PRECISION YV,WW(50) DOUBLE PRECISION SKWDX,SKWDY,SKWOF,SUM DOUBLE PRECISION DEGREE,DD,DSIN,DCOS DOUBLE PRECISION YLOW, YHIGH, WSUM C LOGICAL IF1D C CHARACTER*1 METH C DOUBLE PRECISION IX_R8R8, XI_R8R8 C C DEGREE = RAD / DEGREE C DATA DEGREE/0.017453292519943295769237D+0/ C NWW = 50 C C ... define extraction constants C WIDTH = PAR(1) IWIDTH = WIDTH HWIDTH = 0.5*WIDTH SKEW = PAR(2) SKWDX = 0.0 SKWDY = 1.0 SKWOF = 0.0 XMIN = 0.0 XMAX = XMIN C OFFSET = PAR(3) - HWIDTH + 0.5 OFFSET = PAR(3) IF (SKEW.NE.0.0) THEN DD = SKEW*DEGREE SKWDX = -DSIN(DD) ! DX PER PIXEL SKWDY = DCOS(DD) ! DY PER PIXEL SKWOF = -SKWDX*HWIDTH ! 2D INTERPOLATION IF1D = .FALSE. END IF AX = 1.D0 + SKWOF IYLIM = NPIXA2 - IWIDTH NX = SKWOF N1 = NX N2 = NPIXA1 - NX C C ... get ordinate values C K = PAR(8) L = PAR(9) YV = PAR(10) C C ... loop over the samples C DO 10 IX = 1,N1 OUT(IX,INDX) = 0. 10 CONTINUE DO 100 IX = N2 + 1,NPIXA1 OUT(IX,INDX) = 0. 100 CONTINUE DO 50 IX1 = N1 + 1,N2 IP = 0 DC = 1.D0 YVAL = 0.D0 XVAL = IX_R8R8(AX,STARTA(1),STEPA(1)) C Computes position of the center of the slit from the coefficients of C the bivariate polynomial stored in table order. DO 30 L1 = 0,L IP = IP + 1 WW(IP) = DC YVAL = YVAL + WW(IP)*DPAR(IP) DO 20 K1 = 1,K IP = IP + 1 WW(IP) = WW(IP-1)*XVAL YVAL = YVAL + WW(IP)*DPAR(IP) 20 CONTINUE DC = DC*YV 30 CONTINUE C C Corrects the central position for the offset C OFFSET does not take into account HWIDTH C YVAL = XI_R8R8(YVAL,STARTA(2),STEPA(2)) + OFFSET C YLOW = YVAL - HWIDTH YHIGH = YVAL + HWIDTH LOSLIT = INT(YLOW) HISLIT = INT(YHIGH) IF (LOSLIT.LT.1) LOSLIT = 1 IF (HISLIT.GT.NPIXA2) HISLIT = NPIXA2 SUM = 0.D0 WSUM = 0.D0 IX = INT(AX) C IF (LOSLIT.EQ.HISLIT) THEN SUM = (YHIGH-YLOW)*IN(IX,LOSLIT) C ELSE C DO 40 IW = LOSLIT,HISLIT IF (IW.EQ.LOSLIT) THEN SUM = SUM + IN(IX,IW)*(LOSLIT+1.-YLOW) ELSE IF (IW.EQ.HISLIT) THEN SUM = SUM + IN(IX,IW)*(YHIGH-HISLIT) ELSE SUM = SUM + IN(IX,IW) ENDIF 40 CONTINUE C ENDIF C IF (METH(1:1) .EQ. 'L') THEN ! Linear extraction OUT(IX1,INDX) = SUM ELSE ! Average extraction OUT(IX1,INDX) = SUM/WIDTH ENDIF XMIN = AMIN1(XMIN,OUT(IX1,INDX)) XMAX = AMAX1(XMAX,OUT(IX1,INDX)) AX = AX + 1.D0 50 CONTINUE RETURN END SUBROUTINE EXTR2M(IN,NPIXA1,NPIXA2,OUT,NB1,NB2,COL,NWDTH,INDX,NPAR $ ,PAR,ND,DPAR,XMIN,XMAX,METH,STARTA,STEPA) C C EXTRACTION ROUTINE - MEDIAN C C ARGUMENTS C IN REAL*4 INPUT IMAGE C NPIXA1 INTG*4 NO. OF SAMPLES C NPIXA2 INTG*4 NO. OF LINES C OUT REAL*4 EXTRACTED ORDERS C NB1 INTG*4 NO. OF SAMPLES C NB2 INTG*4 NO. OF ORDERS C INDX INTG*4 SEQUENTIAL NO. OF THE ORDER C NPAR INTG*4 NO. OF EXTRACTION PARAMETERS C PAR REAL*4 EXTRACTION PARAMETERS C ND INTG*4 NO. OF COEFFS. C DPAR REAL*8 COEFFS. TO DEFINE ORDER POSITION C XMIN REAL*4 MINIMUM VALUE C XMAX REAL*4 MAXIMUMVALUE C IMPLICIT NONE C INTEGER NB1, NB2, NPAR, NPIXA1, NPIXA2, NWDTH, NWDTH2, ND INTEGER NWW, IWIDTH, IYLIM, NX, N1, N2, K, L, IX, INDX INTEGER IX1, IP, L1, K1, IW INTEGER LOSLIT, HISLIT C REAL OUT(NB1,NB2),PAR(NPAR) REAL IN(NPIXA1,NPIXA2),COL(NWDTH) REAL XMIN,XMAX,LIMITS(2) C DOUBLE PRECISION DPAR(ND),STARTA(2),STEPA(2) DOUBLE PRECISION XVAL,YVAL,AX DOUBLE PRECISION WIDTH,SKEW,HWIDTH,OFFSET,DC DOUBLE PRECISION YV,WW(50) DOUBLE PRECISION SKWDX,SKWDY,SKWOF,SUM DOUBLE PRECISION DEGREE,DD,DSIN,DCOS DOUBLE PRECISION YLOW, YHIGH, WSUM DOUBLE PRECISION IX_R8R8, XI_R8R8 C LOGICAL IF1D C CHARACTER*1 METH CHARACTER*80 OUTPUT C C DEGREE = RAD / DEGREE C DATA DEGREE/0.017453292519943295769237D+0/ DATA LIMITS/0.0,0.0/ C NWW = 50 C C ... define extraction constants C WIDTH = PAR(1) IWIDTH = WIDTH HWIDTH = 0.5*WIDTH NWDTH2 = HWIDTH SKEW = PAR(2) SKWDX = 0.0 SKWDY = 1.0 SKWOF = 0.0 XMIN = 0.0 XMAX = XMIN C OFFSET = PAR(3) - HWIDTH + 0.5 OFFSET = PAR(3) IF (SKEW.NE.0.0) THEN DD = SKEW*DEGREE SKWDX = -DSIN(DD) ! DX PER PIXEL SKWDY = DCOS(DD) ! DY PER PIXEL SKWOF = -SKWDX*HWIDTH ! 2D INTERPOLATION IF1D = .FALSE. END IF AX = 1.D0 + SKWOF IYLIM = NPIXA2 - IWIDTH NX = SKWOF N1 = NX N2 = NPIXA1 - NX C C ... get ordinate values C K = PAR(8) L = PAR(9) YV = PAR(10) C C ... loop over the samples C DO 10 IX = 1,N1 OUT(IX,INDX) = 0. 10 CONTINUE DO 100 IX = N2 + 1,NPIXA1 OUT(IX,INDX) = 0. 100 CONTINUE DO 50 IX1 = N1 + 1,N2 IP = 0 DC = 1.D0 YVAL = 0.D0 XVAL = IX_R8R8(AX,STARTA(1),STEPA(1)) C Computes position of the center of the slit from the coefficients of C the bivariate polynomial stored in table order. DO 30 L1 = 0,L IP = IP + 1 WW(IP) = DC YVAL = YVAL + WW(IP)*DPAR(IP) DO 20 K1 = 1,K IP = IP + 1 WW(IP) = WW(IP-1)*XVAL YVAL = YVAL + WW(IP)*DPAR(IP) 20 CONTINUE DC = DC*YV 30 CONTINUE C C Corrects the central position for the offset C OFFSET does not take into account HWIDTH C YVAL = XI_R8R8(YVAL,STARTA(2),STEPA(2)) + OFFSET C YLOW = YVAL - HWIDTH YHIGH = YVAL + HWIDTH LOSLIT = INT(YLOW) HISLIT = INT(YHIGH) IF (LOSLIT.LT.1) LOSLIT = 1 IF (HISLIT.GT.NPIXA2) HISLIT = NPIXA2 SUM = 0.D0 WSUM = 0.D0 IX = INT(AX) C IF (LOSLIT.EQ.HISLIT) THEN OUT(IX1,INDX) = IN(IX,LOSLIT) ELSE K1 = 0 DO 40 IW = LOSLIT,HISLIT K1 = K1 + 1 IF (IW.EQ.LOSLIT) THEN COL(K1) = IN(IX,IW)!*(LOSLIT+1.-YLOW) ELSE IF (IW.EQ.HISLIT) THEN COL(K1) = IN(IX,IW)!*MIN(1.0,YHIGH-HISLIT) ELSE COL(K1) = IN(IX,IW) ENDIF 40 CONTINUE !determine MEDIAN: L1 = K1 / 2 IF (K1.GT.0) THEN CALL GET_MEDIAN(COL,NWDTH,K1,L1,OUT(IX1,INDX),IW) IF (IW.LT.0) THEN WRITE(OUTPUT,12345) IX1,INDX CALL STTPUT(OUTPUT) ENDIF ELSE OUT(IX1,INDX) = -1.0 ENDIF ENDIF C XMIN = AMIN1(XMIN,OUT(IX1,INDX)) XMAX = AMAX1(XMAX,OUT(IX1,INDX)) AX = AX + 1.D0 50 CONTINUE C RETURN C 12345 FORMAT('ERROR at ',I5,',',I5,' occurred') C END SUBROUTINE GET_MEDIAN(RA,N,NUSE,LL,VALU,STATUS) IMPLICIT NONE C INTEGER N,LL,STATUS INTEGER NUSE C REAL RA(N),VALU C VALU = 0.0 STATUS = 0 C IF (NUSE.LE.3) THEN IF (NUSE.LT.1) THEN STATUS = -1 RETURN ELSE IF (NUSE.LE.2) THEN STATUS = 2 VALU = RA(1) c VALU = (RA(1) + RA(2))/2.0 !for ndim = 1,2, median is 1st elem. RETURN ELSE STATUS = 3 VALU = RA(2) !for ndim = 3, median is 2nd elem. RETURN ENDIF ENDIF C CALL NSORT(RA,NUSE,LL,VALU) RETURN C END