C @(#)necrebi.for 17.1.1.1 (ESO-IPG) 01/25/02 17:51:31 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 PROGRAM ECHREB C C+ C.COPYRIGHT: Copyright (c) 1987 European Southern Observatory, C all rights reserved C C.VERSION: 1.0 ESO-FORTRAN Conversion, AA 10:37 - 11 JAN 1988 C C.LANGUAGE: F77+ESOext C C.AUTHOR: D.PONZ C C C.IDENTIFICATION C C program ECHREBI1 C C.MODIFICATIONS C C 910128 P.Ballester Define and set variables C C.KEYWORDS C C ECHELLE, CASPEC, REBIN C C.PURPOSE C C OBTAIN A 1D FRM CALIBRATED IN WAVELENGTHS FROM A 2D FRM C OF TYPE ECHELLE WITH THE DISPERSION COEFFS DEFINED AS DESCRIPTORS C C.ALGORITHM C C EACH LINE IN THE INPUT FRM IS REBINNED IN LINEAR STEP IN WAVELENGTH C ACCORDING TO THE DISPERSION COEFFS. C C.INPUT/OUTPUT C C REBIN/ECHELLE INPUT OUTPUT STEP [METHOD] [TABLE] C OR C REBIN/ECHELLE INPUT OUTPUT REFER_IMAGE [METHOD] [TABLE] C C METHOD CAN BE EITHER "NONLINEAR" - NON LINEAR REBINNING USING C THE COEFFS IN THE TABLE C "LINEAR" - LINEAR REBINNING C C.MODIFICATIONS C [2.2] DEFINE THE SAMPLING DOMAIN VIA A REFER IMAGE C C 010201 last modif C C C- IMPLICIT NONE C INTEGER MADRID, ID, IR INTEGER I,I1,IMNO INTEGER*8 IPNTRA INTEGER*8 IPNTRB INTEGER ISTAT,NAXISA,NAXISB INTEGER PIXBN, KUN,KNUL, NNNN INTEGER NPIXA(3),NPIXB(3),IORD(2),IORD2(2),TID INTEGER IORDER(100),NCOE(2),IA1(100),NPTOT(100) DOUBLE PRECISION WSTART(100),STEPA(3),STARTA(3) DOUBLE PRECISION STEPB(3),STARTB(3) DOUBLE PRECISION WSTEP, A1(7,100) REAL CUT(4),XMIN,XMAX CHARACTER*60 OUTFRM,INFRM,REFER,TABLE CHARACTER METHOD*1, LINE*80, DSCNAM*8 CHARACTER IDENTA*72,CUNITA*64,CUNITB*64 LOGICAL IREF INCLUDE 'MID_INCLUDE:ST_DEF.INC' COMMON /VMR/MADRID(1) INCLUDE 'MID_INCLUDE:ST_DAT.INC' DATA CUNITB/'Flux Relative wavel.'/ DATA NPIXA /3*1/, STARTA /3*0.0/, STEPA /3*1.0/ DATA NPIXB /3*1/, STARTB /3*0.0/, STEPB /3*1.0/ C C ... INITIALIZE SYSTEM C CALL STSPRO('ECHREB') CALL STKRDC('IN_A',1,1,60,I,INFRM,KUN,KNUL,ISTAT) CALL STKRDC('OUT_A',1,1,60,I,OUTFRM,KUN,KNUL,ISTAT) CALL CLNFRA(INFRM,INFRM,0) CALL CLNFRA(OUTFRM,OUTFRM,0) CALL STKRDC('P3',1,1,60,I,REFER,KUN,KNUL,ISTAT) CALL STKRDC('P4',1,1,1,I,METHOD,KUN,KNUL,ISTAT) I1 = INDEX('.+1234567890',REFER(1:1)) IF (I1.NE.0) THEN CALL GENCNV(REFER,4,1,NNNN,XMIN,WSTEP,ISTAT) IF (ISTAT.LT.1) CALL STETER(11,'Wrong sampling step') C WRITE(LINE,9010) WSTEP CALL STTPUT(LINE,ISTAT) IREF = .FALSE. ELSE IREF = .TRUE. CALL CLNFRA(REFER,REFER,0) END IF IF (METHOD.EQ.'N' .OR. METHOD.EQ.'n') THEN CALL STKRDC('IN_B',1,1,60,I,TABLE,KUN,KNUL,ISTAT) CALL CLNTAB(TABLE,TABLE,0) CALL TBTOPN(TABLE,F_I_MODE,TID,ISTAT) CALL STDRDI(TID,'ORDER',1,2,I,IORD,KUN,KNUL,ISTAT) CALL STDRDI(TID,'COEFS',1,2,I,NCOE,KUN,KNUL,ISTAT) CALL STDRDD(TID,'COEFD',1,7*NCOE(2),I,A1, . KUN,KNUL,ISTAT) CALL STDRDI(TID,'COEFI',1,NCOE(2),I,IA1, . KUN,KNUL,ISTAT) CALL TBTCLO(TID,ISTAT) END IF C C ... MAP INPUT IMAGE C CALL STIGET(INFRM,D_R4_FORMAT,F_I_MODE,F_IMA_TYPE, . 3,NAXISA,NPIXA,STARTA,STEPA,IDENTA,CUNITA, . IPNTRA,IMNO,ISTAT) IF (NAXISA.EQ.1 .OR. NPIXA(2).EQ.1) THEN NPIXA(2) = 1 NAXISA = 1 END IF IF (NPIXA(2).GT.100) THEN CALL STETER(9, + 'Module ECHRE31: Too many orders (limit 100). Buffer overflow') ENDIF IF (STEPA(2).NE.1.) THEN CALL STETER(9,'Module ECHREBI: Wrong step(2) .ne. 1 ') ENDIF IORD2(1) = STARTA(2) ! Starting relative order number IORD2(2) = STARTA(2) + NPIXA(2) - 1 ! Final relative order number C C ... MAP OUTPUT FRM C IF (METHOD.EQ.'N' .OR. METHOD.EQ.'n') THEN NAXISB = 2 IF (IREF) THEN CALL STFOPN(REFER,D_R4_FORMAT,0,F_IMA_TYPE, . IR,ISTAT) CALL STDRDI(IR,'NPIX',1,2,I,NPIXB,KUN,KNUL,ISTAT) CALL STDRDD(IR,'START',1,2,I,STARTB,KUN,KNUL,ISTAT) CALL STDRDD(IR,'STEP',1,2,I,STEPB,KUN,KNUL,ISTAT) WSTEP = STEPB(1) CALL STDRDI(IR,'NORDER',1,NPIXB(2),I,IORDER, . KUN,KNUL,ISTAT) CALL STDRDI(IR,'NPTOT',1,NPIXB(2),I,NPTOT, . KUN,KNUL,ISTAT) CALL STDRDD(IR,'WSTART',1,NPIXB(2),I,WSTART, . KUN,KNUL,ISTAT) DO 10 I = 1,NPIXB(2) WSTART(I) = WSTART(I)/1000. 10 CONTINUE ELSE CALL WRANGE(NCOE(2),A1,IA1,WSTEP,IORD,NPIXA(1),WSTART, . NPIXB(1),NPTOT,IORDER,STARTA,STEPA,IORD2) NPIXB(2) = NPIXA(2) STARTB(1) = 0.D0 STARTB(2) = STARTA(2) STEPB(1) = WSTEP STEPB(2) = STEPA(2) END IF CALL STIPUT(OUTFRM,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE, . NAXISB,NPIXB,STARTB,STEPB, . IDENTA,CUNITB,IPNTRB,ID,ISTAT) C C ... NON LINEAR REBINNING C CALL ECHRE3(MADRID(IPNTRA),NPIXA(1),NPIXA(2),STARTA,STEPA, . NCOE(2),A1,IA1,IORD,MADRID(IPNTRB),NPIXB(1), . NPIXB(2),WSTART,WSTEP,IORD2) C C ... WRITE PROCESS DESCRIPTORS C NNNN = NPIXB(1)*NPIXB(2) CALL MNMX(MADRID(IPNTRB),NNNN,XMIN,XMAX) CUT(1) = XMIN CUT(2) = XMAX CUT(3) = XMIN CUT(4) = XMAX DO 20 I = 1,NPIXB(2) WSTART(I) = 1000.*WSTART(I) 20 CONTINUE CALL STDCOP(IMNO,ID,3,DSCNAM,ISTAT) CALL STDWRD(ID,'WSTART',WSTART,1,NPIXB(2),KUN,ISTAT) CALL STDWRI(ID,'NPTOT',NPTOT,1,NPIXB(2),KUN,ISTAT) CALL STDWRI(ID,'NORDER',IORDER,1,NPIXB(2),KUN,ISTAT) CALL STDWRR(ID,'LHCUTS',CUT,1,4,KUN,ISTAT) ELSE C ... LINEAR REBINNING IF (NPIXA(2).GT.1) THEN CCCCCCCCCCCCC CALL STDFND(IMNO,'WSTART',I1,I2,I3,ISTAT) C IF (TYPE(1:1).EQ.' ') THEN C CALL STTPUT(' Wrong input image ... ',ISTAT) C GO TO 40 CCCCCCCCCCCCC END IF CALL STDRDD(IMNO,'WSTART',1,NPIXA(2),I,WSTART, . KUN,KNUL,ISTAT) CALL STDRDI(IMNO,'NPTOT',1,NPIXA(2),I,NPTOT, . KUN,KNUL,ISTAT) END IF CALL STDRDR(IMNO,'LHCUTS',1,4,I,CUT,KUN,KNUL,ISTAT) PIXBN = STEPA(1)/WSTEP NAXISB = NAXISA NPIXB(1) = NPIXA(1)*PIXBN NPIXB(2) = NPIXA(2) STARTB(1) = STARTA(1) STARTB(2) = STARTA(2) STEPB(1) = WSTEP STEPB(2) = STEPA(2) IF (NPIXB(2).GT.1) THEN DO 30 I = 1,NPIXB(2) NPTOT(I) = NPTOT(I)*PIXBN 30 CONTINUE ELSE WSTART(1) = STARTB(1) NPTOT(1) = NPIXB(1) END IF CALL STIPUT(OUTFRM,D_R4_FORMAT,F_O_MODE,F_IMA_TYPE, . NAXISB,NPIXB,STARTB,STEPB,IDENTA, . CUNITB,IPNTRB,ID,ISTAT) CALL ECHR01(MADRID(IPNTRA),NPIXA(1),NPIXA(2),STARTA,STEPA, . NPTOT,WSTART, . MADRID(IPNTRB),NPIXB(1),NPIXB(2),WSTEP) CALL STDCOP(IMNO,ID,3,DSCNAM,ISTAT) IF (NPIXB(2).GT.1) THEN CALL STDWRD(ID,'WSTART',WSTART,1,NPIXB(2),KUN,ISTAT) CALL STDWRI(ID,'NPTOT',NPTOT,1,NPIXB(2),KUN,ISTAT) END IF CALL STDWRR(ID,'LHCUTS',CUT,1,4,KUN,ISTAT) END IF CALL STSEPI C10000 FORMAT(G.0) 9000 FORMAT (G6.0) 9010 FORMAT(' Sampling step = ',F10.4) END SUBROUTINE ECHRE3(X,NPIXA1,NPIXA2,STARTA,STEPA,NO,A,IA,IORD, + Y,NPIXB1,NPIXB2,YSTR,YSTP,IORD2) C C REBIN THE ORDERS C IMPLICIT NONE C INTEGER NO, NPIXB1, NPIXB2 INTEGER NPIXA1,NPIXA2,IORD(2),IA(NO),IORD2(2) REAL X(NPIXA1,NPIXA2) REAL Y(NPIXB1,NPIXB2) DOUBLE PRECISION YSTR(1), YSTP, STARTA(2), STEPA(2) DOUBLE PRECISION A(7,NO),B1(3) C INTEGER KMIN, KMAX, I0, I2, I1, K1, II, MIN, MAX, I3 DOUBLE PRECISION WLSTP, WLST2 C C ... GET LIMITS FOR THE FIRST ORDER (LOWER WAVELENGTHS) C KMIN = MIN(IORD(1),IORD(2)) KMAX = MAX(IORD(1),IORD(2)) IF (IORD(1).LT.IORD(2)) THEN I0 = NO I2 = -1 ELSE I0 = 1 I2 = 1 END IF I1 = I0 WLSTP = YSTP/1000. C C ... ITERATION ON ORDERS C DO 20 K1 = IORD2(2), IORD2(1), -1 I3 = I1 + IORD2(1) - 1 IF (IA(I3).LT.0) THEN B1(1) = A(1,I3) B1(2) = -A(2,I3) B1(3) = A(3,I3) END IF WLST2 = YSTR(I1) IF (IA(I3).LT.0) THEN CALL LREBIN(X(1,I1),NPIXA1,STARTA,STEPA,B1,3,Y(1,I1), + NPIXB1,WLST2,WLSTP) C C ... SIGN CHANGE DO 10 II = 1,NPIXB1 Y(II,I1) = -Y(II,I1) 10 CONTINUE ELSE CALL LREBI1(X(1,I1),NPIXA1,STARTA,STEPA,A(1,I3),IA(I3), + Y(1,I1),NPIXB1,WLST2,WLSTP) END IF I1 = I1 + I2 IF (I1.GT.NO) RETURN 20 CONTINUE RETURN END DOUBLE PRECISION FUNCTION POLEV(IA,A,XX) C C C IMPLICIT NONE C INTEGER IA DOUBLE PRECISION A(IA),XX,YY YY = A(1) + XX* (A(2)+XX* (A(3)+XX* (A(4)+XX* (A(5)+XX* (A(6)+ + XX*A(7)))))) POLEV = YY RETURN END SUBROUTINE MNMX(X,N,XMIN,XMAX) C C FIND MIN MAX VALS C IMPLICIT NONE INTEGER N REAL X(N), XMIN, XMAX C INTEGER I C XMIN = X(1) XMAX = XMIN DO 10 I = 2,N XMIN = AMIN1(XMIN,X(I)) XMAX = AMAX1(XMAX,X(I)) 10 CONTINUE RETURN END SUBROUTINE WRANGE(NCOE,A,IA,WSTEP,IORD,NPIXA,WSTART,NPIXB,NPTOT, + IORDER,STARTA,STEPA,IORD2) C C GET STARTING POSITIONS FOR EACH ORDER AND NO. OF PIXELS PER ORDER C C C ... GET LIMITS FOR THE FIRST ORDER (LOWER WAVELENGTHS) C IMPLICIT NONE C INTEGER NCOE, NPIXA, NPIXB DOUBLE PRECISION A(7,NCOE),POLEV,B(3), STARTA(1), STEPA(1) DOUBLE PRECISION WSTART(1), WSTEP, WLST1, WLEN1 INTEGER IA(NCOE),IORD(2),NPTOT(1),IORDER(1), IORD2(2) C INTEGER KMIN, KMAX, I0, I1, I2, I, ISTEP, K, N1, I3 c INTEGER*8 N1 DOUBLE PRECISION WS1, WW1, IX_IR8 C KMIN = MIN(IORD(1),IORD(2)) KMAX = MAX(IORD(1),IORD(2)) IF (IORD(1).LT.IORD(2)) THEN I0 = NCOE I2 = -1 ELSE I0 = 1 I2 = 1 END IF I1 = I0 NPIXB = 0 WS1 = WSTEP/1000. C C ... get order number C ISTEP = 1 IF (IORD(1).GT.IORD(2)) ISTEP = -1 ! Default case in present versions IORDER(1) = IORD(1) - IORD2(1) + 1 DO 10 I = 2, IORD2(2) - IORD2(1) + 1 IORDER(I) = IORDER(I-1) + ISTEP 10 CONTINUE C C ... iteration on orders to get starting wavelength and no. of pixels C DO 20 K = IORD2(2), IORD2(1), -1 I3 = I1 + IORD2(1) - 1 IF (IA(I3).LT.0) THEN B(1) = A(1,I3) B(2) = -A(2,I3) B(3) = A(3,I3) CALL ECHORD(B,3,NPIXA,WLST1,WLEN1,STARTA,STEPA) ELSE WW1 = POLEV(IA(I3),A(1,I3),IX_IR8(1,STARTA,STEPA)) WLST1 = WW1 WW1 = POLEV(IA(I3),A(1,I3),IX_IR8(NPIXA,STARTA,STEPA)) WLEN1 = WW1 END IF WLST1 = WLST1*10. WLEN1 = WLEN1*10. N1 = (WLEN1-WLST1)/WS1 IF (N1.GT.65534) THEN !added swolf@eso.org N1 = 0 !if WLEN1 and WLST1 are NaN N1 gets wrong! ENDIF !--> NPIXB gets wrong, too! NPTOT(I1) = N1 NPIXB = MAX(NPIXB,N1) IF (I1.EQ.I0) THEN WSTART(I1) = WS1*NINT(WLST1/WS1) !before: WLST1 (swolf@eso.org) ELSE WSTART(I1) = WSTART(I0) + WS1*NINT((WLST1-WSTART(I0))/WS1) END IF I1 = I1 + I2 IF (I1.GT.NCOE) GO TO 30 20 CONTINUE 30 NPIXB = NPIXB + 3 RETURN END SUBROUTINE ECHR01(A,NPIXA1,NPIXA2,STARTA,STEPA,NPTOT,WSTART, + B,NPIXB1,NPIXB2,WSTEP) C C REBIN ECHELLE ORDERS C IMPLICIT NONE C INTEGER NPIXA1,NPIXA2,NPIXB1,NPIXB2 DOUBLE PRECISION STEPA(2),STARTA(2) REAL A(NPIXA1,NPIXA2),B(NPIXB1,NPIXB2) DOUBLE PRECISION WSTART(NPIXA2), WSTEP INTEGER NPTOT(NPIXA2) C INTEGER NORD, IORD, I, NPREB DOUBLE PRECISION XSTART, XP, ADEX, PIXBN REAL VBIN C PIXBN = WSTEP/STEPA(1) NORD = NPIXA2 DO 30 IORD = 1,NORD XSTART = WSTART(IORD) NPREB = NPTOT(IORD) DO 10 I = 1,NPREB XP = XSTART + (I-1)*WSTEP ADEX = (XP-XSTART)/STEPA(1) + 1. B(I,IORD) = VBIN(A(1,IORD),PIXBN,ADEX) 10 CONTINUE DO 20 I = NPREB + 1,NPIXB1 B(I,IORD) = 0. 20 CONTINUE 30 CONTINUE RETURN END REAL FUNCTION VBIN(ARRAY,PIXPBN,ADEX) C C FUNCTION TO RETURN A VALUE FOR A REBINNED BIN C C INPUT ARGUMENTS C ARRAY REAL*4 ARRAY WITH REAL ORIGINAL PIXELS C PIXPBN REAL*4 PIXELS PER BIN, NUMBER OF ORIGINAL C PIXELS IN THE NEW BIN, WITH FRACTIONS C ADEX REAL*4 CENTER INDEX IN THE ARRAY FOR THE NEW BIN, C WITH FRACTIONS C IMPLICIT NONE INTEGER IBOTDX,ITOPDX,IENLP,I,ISTLP REAL HBIN,BOTDX,TOPDX,SUM,BOTFRA,TOPFRA REAL ARRAY(1) DOUBLE PRECISION PIXPBN, ADEX C ! HALF BIN HBIN = PIXPBN*0.5 ! BOTTOM POINTER, FRACTION BOTDX = ADEX - HBIN ! TOP POINTER, FRACTION TOPDX = ADEX + HBIN ! BOTTOM POINTER, WITHOUT FRACTION IBOTDX = BOTDX ! TOP POINTER, WITHOUT FRACTION ITOPDX = TOPDX C IF (IBOTDX.EQ.ITOPDX) THEN C C BIN SMALLER THAN ORIGINAL PIXELS C AND NEW BIN WITHIN AN ORIGINAL PIXEL C SUM = PIXPBN*ARRAY(IBOTDX+1) ELSE C C BIN LARGER THAN ORIGINAL PIXELS, OR C SMALL BIN OVER TWO PIXELS C BOTFRA = 1. - (BOTDX-FLOAT(IBOTDX)) TOPFRA = TOPDX - FLOAT(ITOPDX) ! START LOOP INDEX ISTLP = IBOTDX + 2 IENLP = ITOPDX IBOTDX = IBOTDX + 1 ITOPDX = ITOPDX + 1 C SUM = 0. IF (IENLP.GE.ISTLP) THEN DO 10 I = ISTLP,IENLP SUM = SUM + ARRAY(I) 10 CONTINUE END IF ! LOWER FRACTION SUM = SUM + BOTFRA*ARRAY(IBOTDX) ! UPPER FRACTION SUM = SUM + TOPFRA*ARRAY(ITOPDX) END IF C C AREA CORRECTION C VBIN = SUM/PIXPBN RETURN END SUBROUTINE ECHORD(AC,NC,NPIX,W1,W2,STARTA,STEPA) C C C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C.COPYRIGHT: Copyright (c) 1987 European Southern Observatory, C all rights reserved C C.VERSION: 1.0 ESO-FORTRAN Conversion, AA 22:34 - 3 DEC 1987 C C.LANGUAGE: F77+ESOext C C.AUTHOR: D.PONZ C C.KEYWORDS: C SPECTROSCOPY, CASPEC, ORDER DEFINITION C C.PURPOSE: C GET WAVELENGTH RANGE OF A GIVEN ORDER C C.ALGORITHM: C C.INPUT/OUTPUT C INPUT ARGUMENTS C C AC(NC) REAL*8 COEFFICIENTS OF THE TRANSFORMATION C NC REAL*4 NO. OF COEFFICIENTS C NPIX INTG*4 NO OF SAMPLES C C OUTPUT ARGUMENTS C C W1 REAL*4 LOWER LIMIT C W2 REAL*4 UPPER LIMIT C C------------------------------------------------------------- C C IMPLICIT NONE INTEGER NPIX, NC DOUBLE PRECISION AC(NC), STARTA(1), STEPA(1) DOUBLE PRECISION W1, W2, IX_IR8 C DOUBLE PRECISION A0,A1,A2,F1,FN C C C ... ASSIGN COEFFS C C F1 = 1.D0 C FN = NPIX F1 = IX_IR8(1,STARTA(1),STEPA(1)) FN = IX_IR8(NPIX,STARTA(1),STEPA(1)) A0 = AC(1) A1 = AC(2) IF (NC.GE.3) THEN A2 = AC(3) W1 = (-A1+DSQRT(A1*A1-4.D0*A2* (A0-F1)))/ (2.D0*A2) W2 = (-A1+DSQRT(A1*A1-4.D0*A2* (A0-FN)))/ (2.D0*A2) ELSE W1 = (F1-A0)/A1 W2 = (FN-A0)/A1 END IF RETURN END SUBROUTINE LREBI1(X,NX,XSTR,XSTP,AC,NC,Y,NY,YSTR,YSTP) C C C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C.COPYRIGHT: Copyright (c) 1987 European Southern Observatory, C all rights reserved C C.VERSION: 1.0 ESO-FORTRAN Conversion, AA 0:22 - 4 DEC 1987 C C.LANGUAGE: F77+ESOext C C.AUTHOR: D.PONZ C C.IDENTIFICATION C SUBROUTINE LREBIN1 VERSION 1 22 MAR 1983 C J.D.PONZ C C.KEYWORDS: C SPECTROSCOPY, CASPEC, REBINNING C C.PURPOSE: C DATA REBINNING C CONVERSION FROM SAMPLES INTO WAVELENGTHS C C.ALGORITHM: C LINEAR INTERPOLATION. THE TRANSFORMATION SAMPLE NUMBER INTO C WAVELENGTH IS DEFINED BY C W = A0 + A1*S + A2*S**2 + ... C C.INPUT/OUTPUT C INPUT ARGUMENTS C C X(NX) REAL*4 INPUT ARRAY C NX INTG*4 NO OF SAMPLES C XSTR REAL*4 STARTING VALUE OF INPUT C XSTP REAL*4 INPUT STEP C AC(NC) REAL*8 COEFFICIENTS OF THE TRANSFORMATION C NC REAL*4 NO. OF COEFFICIENTS C NY INTG*4 NO. OF OUTPUT VALUES C YSTR REAL*4 STARTING VALUE OF THE INPUT C YSTP REAL*4 STEP OF THE OUTPUT C C OUTPUT ARGUMENTS C C Y(NY) REAL*4 OUTPUT ARRAY C C------------------------------------------------------------- C C IMPLICIT NONE INTEGER NX, NY, NC REAL X(NX),Y(NY) DOUBLE PRECISION YSTR, YSTP, XSTR, XSTP C DOUBLE PRECISION AC(NC),A0,A1,A2,A3,A4,A5,A6,FN,FL,WPOLY DOUBLE PRECISION FACTOR,DPOLY,CNORM, IX_R8R8, XN, XNP1 DOUBLE PRECISION PN,PNP1,TM,TMP1,C,COV,VAL, COVX DOUBLE PRECISION B2, B3, B4, B5, B6, DELT, STARTL INTEGER M1, MST, N, K, MPT C WPOLY(FL) = A0 + FL* (A1+FL* (A2+FL* (A3+FL* (A4+FL* (A5+FL* + A6))))) DPOLY(FL) = A1 + FL* (B2+FL* (B3+FL* (B4+FL* (B5+FL*B6)))) C C ... ASSIGN COEFFS C A0 = AC(1) A1 = AC(2) A2 = AC(3) A3 = AC(4) A4 = AC(5) A5 = AC(6) A6 = AC(7) B2 = 2.D0*A2 B3 = 3.D0*A3 B4 = 4.D0*A4 B5 = 5.D0*A5 B6 = 6.D0*A6 CNORM = 10.D0/YSTP C C ... CHECK FOR OUTPUT START C STARTL = YSTR - 0.5*YSTP FN = 0.5D0 XN = IX_R8R8(FN,XSTR,XSTP) PN = 10.D0*WPOLY(XN) DELT = PN - STARTL IF (DELT.GT.0.) THEN M1 = DELT/YSTP + 1. TMP1 = PN + YSTP TM = PN ELSE M1 = 1 TMP1 = STARTL + YSTP TM = STARTL END IF C C ... CHECK FOR INPUT START C DO 10 N = 1,NX FN = FN + 1.D0 XNP1 = IX_R8R8(FN,XSTR,XSTP) PNP1 = 10.D0*WPOLY(XNP1) IF (PNP1.GT.TM) GO TO 20 PN = PNP1 10 CONTINUE N = NX + 1 GO TO 30 20 C = X(N) 30 VAL = 0. COVX = XNP1 - XN COV = (PNP1 - PN)/COVX C C ... REBINNING LOOP C MPT = 1 MST = 1 K = 0 35 CONTINUE IF ( .NOT. (MPT.LE.NY)) GO TO 80 K = K + 1 IF ((K.LT.M1) .OR. (N.GT.NX)) GO TO 60 C 40 IF (PNP1.EQ.PN) GO TO 90 VAL = VAL + C* (DMIN1(TMP1,PNP1)-DMAX1(TM,PN))/COV IF (TMP1.LT.PNP1) GO TO 50 N = N + 1 IF (N.GT.NX) GO TO 60 PN = PNP1 XN = XNP1 FN = FN + 1.D0 XNP1 = IX_R8R8(FN,XSTR,XSTP) PNP1 = 10.D0*WPOLY(XNP1) C = X(N) COV = (PNP1 - PN)/(XNP1-XN) GO TO 40 50 TM = TMP1 TMP1 = TMP1 + YSTP 60 FACTOR = CNORM*DPOLY(IX_R8R8(FN,XSTR,XSTP)) C FACTOR = FACTOR*10.D0*DABS(DPOLY(IX_R8R8(FN,XSTR,XSTP)))/YSTP Y(MPT) = VAL*FACTOR MPT = MPT + 1 VAL = 0. GOTO 35 80 CONTINUE 90 CONTINUE DO 100 K = MPT + 1,NY Y(K) = 0. 100 CONTINUE RETURN END SUBROUTINE LREBIN(X,NX,XSTR,XSTP,AC,NC,Y,NY,YSTR,YSTP) C C C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C.COPYRIGHT: Copyright (c) 1987 European Southern Observatory, C all rights reserved C C.VERSION: 1.0 ESO-FORTRAN Conversion, AA 19:08 - 6 DEC 1987 C C.LANGUAGE: F77+ESOext C C.AUTHOR: D.PONZ C C.KEYWORDS: C SPECTROSCOPY, CASPEC, REBINNING C C.PURPOSE: C DATA REBINNING C CONVERSION FROM SAMPLES INTO WAVELENGTHS C C.ALGORITHM: C LINEAR INTERPOLATION. THE TRANSFORMATION WAVELENGTH INTO C SAMPLE NUMBER IS DEFINED BY C S = A0+A1*T+A2*T**2 C C.INPUT/OUTPUT C INPUT ARGUMENTS C C X(NX) REAL*4 INPUT ARRAY C NX INTG*4 NO OF SAMPLES C XSTR REAL*4 STARTING VALUE OF INPUT C XSTP REAL*4 INPUT STEP C AC(NC) REAL*8 COEFFICIENTS OF THE TRANSFORMATION C NC REAL*4 NO. OF COEFFICIENTS C NY INTG*4 NO. OF OUTPUT VALUES C YSTR REAL*4 STARTING VALUE OF THE INPUT C YSTP REAL*4 STEP OF THE OUTPUT C C OUTPUT ARGUMENTS C C Y(NY) REAL*4 OUTPUT ARRAY C C------------------------------------------------------------- C C IMPLICIT NONE INTEGER NX, NY, NC, N, M1 REAL X(NX),Y(NY) DOUBLE PRECISION AC(NC),A0,A1,A2,A3,A4,A5,FN DOUBLE PRECISION FACTOR,CNORM, START, PN, DELT, STARTL DOUBLE PRECISION TM, TMP1 DOUBLE PRECISION IX_R8R8 DOUBLE PRECISION XSTR, XSTP, YSTR, YSTP C INTEGER K, MPT, MST DOUBLE PRECISION PNP1, C, VAL, COV, DSQRT, DMIN1, DMAX1 DOUBLE PRECISION XN, XNP1, COVX C C WPOLY(FL) = FL*(A1+FL*(A2+FL*(A3+FL*(A4+FL*A5)))) C C ... ASSIGN COEFFS C A0 = AC(1) A1 = -AC(2) A2 = AC(3) A3 = A1*A1 A4 = 4.D0*A2 A5 = 2.D0*A2 CNORM = 10.D0/ (A5*YSTP) C A0 = AC(1) C A1 = AC(2) C A2 = AC(3) C C ... CHECK FOR OUTPUT START C START = YSTR - 0.5*YSTP C STARTL= START - A0 STARTL = START M1 = 1 FN = 0.5D0 XN = IX_R8R8(FN,XSTR,XSTP) C PN = WPOLY(IX_R8R8(FN,XSTR,XSTP)) PN = 10.D0* (A1+DSQRT(A3-A4* (A0-XN)))/A5 DELT = PN - STARTL IF (DELT.GT.0.) THEN M1 = DELT/YSTP + 1. TMP1 = STARTL + M1*YSTP TM = TMP1 - YSTP ELSE TMP1 = STARTL + YSTP TM = STARTL END IF C C ... CHECK FOR INPUT START C DO 10 N = 1,NX FN = FN + 1.D0 XNP1 = IX_R8R8(FN,XSTR,XSTP) C PNP1 = WPOLY(IX_R8R8(FN,XSTR,XSTP)) PNP1 = 10.D0* (A1+DSQRT(A3-A4* (A0-XNP1)))/A5 IF (PNP1.GT.TM) GO TO 20 PN = PNP1 XN = XNP1 10 CONTINUE N = NX + 1 GO TO 30 20 C = X(N) 30 VAL = 0. COVX = XNP1 - XN COV = (PNP1 - PN)/COVX C C ... REBINNING LOOP C MPT = 1 MST = 1 K = 0 C 35 CONTINUE IF ( .NOT. (MPT.LE.NY)) GO TO 80 K = K + 1 IF ((K.LT.M1) .OR. (N.GT.NX)) GO TO 60 40 VAL = VAL + C* (DMIN1(TMP1,PNP1)-DMAX1(TM,PN))/COV IF (TMP1.LT.PNP1) GO TO 50 N = N + 1 IF (N.GT.NX) GO TO 60 PN = PNP1 XN = XNP1 FN = FN + 1.D0 XNP1 = IX_R8R8(FN,XSTR,XSTP) C PNP1 = WPOLY(IX_R8R8(FN,XSTR,XSTP)) PNP1 = 10.D0* (A1+DSQRT(A3-A4* (A0-XNP1)))/A5 C = X(N) COVX = XNP1 - XN COV = (PNP1 - PN)/COVX GO TO 40 50 TM = TMP1 TMP1 = TMP1 + YSTP 60 FACTOR = CNORM*(DSQRT(A3-A4* (A0-XNP1))- + DSQRT(A3-A4* (A0-XNP1-1.D0))) Y(MPT) = VAL*FACTOR MPT = MPT + 1 VAL = 0. 80 CONTINUE DO 90 K = MPT + 1,NY Y(K) = 0. 90 CONTINUE RETURN END