C @(#)tdrebin.for 17.1.1.1 (ESO-DMD) 01/25/02 17:47:16 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 18:08 - 11 DEC 1987 C C.LANGUAGE: F77+ESOext C C.IDENTIFICATION C TDREBIN.FOR C C C.AUTHOR: M.ROSA C REBPOL modified 900117 M Peron C remove call to NAG 900903 M Peron C C 010222 last modif C C-------------------------------------------------------------------- SUBROUTINE REBMET(NIN,XI,YI,WI,NOUT,XO,WO,IRFUN,NPAR,NPACT,DPAR, + IOPT,MINT,YO,RMIN,RMAX) C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C.KEYWORDS C NonLINEAR rebin C C.PURPOSE C NonLINEAR rebin C C.ALGORITHM C C TBF C C.USAGE C CALL REBMET(NIN,XI,YI,WI,NOUT,XO,WO,IRFUN,NPAR,NPACT,DPAR, C IOPT,MINT,YO,RMIN,RMAX) C C.PARAMETERS C Inputs: C NIN,XI,YI,WI no. of (input) points, indep, dep and bin inputs C NOUT,XO,WO no, of (output) points, indep, bin outputs C IRFUN,NPAR,NPACT,PAR func. no., size par., act no par., parameters C C Note, function relates XO --> XI i.e. xi(f(xo)) !!!!!! C C IRFUN = 1: linear x (old) = A + Bx C = 2: poly. x = A + Bx + Cxx + (<= 16 coeffs) C = 3: inv. x = A + B/x C = 4: exp. x = A + B*EXP(Cx) C = 5: dexp. x = A + B*10.**(Cx) C = 6: log x = A + B*LOG(Cx) C = 7: dlog x = A + B*LOG10(Cx) C = 8: inv.pol. x = A + B/(x-x0) + C/(x-x0)**2 + ... C = 9: user defined function C C IOPT,MINT interpol/integrate option, no. of C points to be interpolated C IOPT = 1: spline interpolation and gaussian integration (most C accurate but most time consuming option) C = 2: intensity is simply relevant portion of pixel's C intensity; for heavy oversampling of input, use this. C = 3: linear interpolation, followed by staightforward C integration (use if i/p and o/p binwidths similar). C C Outputs: C YO,RMIN,RMAX dep and maxmin of output C C.RESTRICTIONS C C.MODIFICATIONS C C.LONG DESCRIPTION C The subroutines below are well commented. For background C information on the use of this routine, see the Users Guide C for the modelling system in USG.MEM in [STMODEL.DOC] C (available through M. Rosa, ST-ECF). C C-------------------------------------------------------------------- C INTEGER NIN,NOUT,IRFUN,NPAR,NPACT,IOPT,MINT INTEGER NA,NZ,I,J,LOP,LIP,IN,KO,KF,LL,LH,ML,MH C DOUBLE PRECISION XI(NIN),YI(NIN),WI(NIN) DOUBLE PRECISION XO(NOUT),WO(NOUT),DPAR(NPAR),PL,PH DOUBLE PRECISION XL,XH,XW,XINCR,XVLS,DLOP,DLIP,DL,DH C REAL YO(NOUT),RMIN,RMAX,VAL,VALI C C C C No. of intervals XVLS = DBLE(MINT-1) C C We start here doing some brain work in order to speed up C the program --- all this should be ok but who knows (MR/19/4/85) C C There are 4 cases for the indep. arrays C XI(i) < XI(i+1) C XI(i) > XI(i+1) C XO(i) < XO(i+1) C XO(i) > XO(i+1) C C and two for the mapping of XO onto XI C XO small onto XI small C XO small onto XI large C C lets make a few tests, so we can optimize the DO LOOPS C XO increasing IF (XO(1).LE.XO(2)) THEN KO = 1 XL = XO(1) XH = XO(NOUT) ELSE KO = -1 XH = XO(1) XL = XO(NOUT) END IF CALL REBFNC(IRFUN,XL,NPAR,NPACT,DPAR) C the function is inverting direction CALL REBFNC(IRFUN,XH,NPAR,NPACT,DPAR) C 1: normal, -1: inverted direction KF = 1 IF (XL.GT.XH) KF = -1 LOP = KF*KO C C So we set up the loops in such a way that the input indep. will be C increasing for all cases. That means we can keep track of search C limits (large gain in search for long arrays). C C REBIN with flux conservation means that we have to integrate the flux C in the input data between the projected flanks of our new pixel. C Remember, that the binwidth is negative if the arrays are C monotonicaly decreasing !!! Since we will have to look for the C projection of the flanks of each pixel XO in the XI array (XL,XR), C where Left and Right referes to the direction of indexing in XI, C we have to be careful with the directions (4 cases): C C XO incr, LOP pos, WO pos, xl = xo - 0.5*wo , xr = xl + wo C XO incr, LOP neg, WO pos, xl = xo + 0.5*wo , xr = xl - wo C XO decr, LOP pos, WO neg, xl = xo + 0.5*wo , xr = xl - wo C XO decr, LOP neg, WO neg, xl = xo - 0.5*wo , xr = xl + wo C C Hence, XL = XO - 0.5*LOP*WO, XR = XL + LOP*WO (in REAL*8) C C The new indep. may not map completely onto the old and vice versa. C We will NOT extrapolate and we will fill undefined output with C NULLVALUE. Now we need C to obtain the applicable limits in XO, the direction in XI and the C search boundaries for XI and fill YO with 0 outside, in order to C avoid all sorts of limit testing in the actual run through XO. C XI increasing ! decreasing LIP = 1 IF (XI(1).GT.XI(2)) LIP = -1 DLIP = 0.5D0*DBLE(LIP) DLOP = 0.5D0*DBLE(LOP) C C Check for the active area in XO C IF (LOP.GT.0) THEN NA = 1 NZ = NOUT ELSE NA = NOUT NZ = 1 ENDIF C DL = XI(1) DH = XI(NIN) C XL = XO(NA) - DLOP*WO(NA) XH = XO(NZ) + DLOP*WO(NZ) CALL REBFNC(IRFUN,XL,NPAR,NPACT,DPAR) CALL REBFNC(IRFUN,XH,NPAR,NPACT,DPAR) C complete mapping of XO IF (XL.GE.DL .AND. XH.LE.DH) GO TO 50 C no mapping of XO IF (XL.GE.DH .OR. XH.LE.DL) THEN WRITE (*,*) ' No overlap between input', + ' and output independent !','/' END IF C check lower DX bound IN = NA DO 10, I=NA,NZ,LOP XL = XO(I) - DLOP*WO(I) IN = I CALL REBFNC(IRFUN,XL,NPAR,NPACT,DPAR) IF (XL.GE.DL) GO TO 20 YO(I) = 0.00 10 CONTINUE 20 NA = IN C check upper DX bound IN = NZ DO 30, I=NZ,NA,-LOP XH = XO(I) + DLOP*WO(I) IN = I CALL REBFNC(IRFUN,XH,NPAR,NPACT,DPAR) IF (XH.LT.DH) GO TO 40 YO(I) = 0.00 30 CONTINUE 40 NZ = IN C the initial lower limit DX 50 IF (LIP.EQ.1) THEN LL = 1 LH = NIN ELSE LL = NIN LH = 1 ENDIF C C ... and finally we go through the XO array C DO 130, I = NA,NZ,LOP VAL = 0.0 XW = DLOP*WO(I) XL = XO(I) - XW XH = XL + XW + XW CALL REBFNC(IRFUN,XL,NPAR,NPACT,DPAR) CALL REBFNC(IRFUN,XH,NPAR,NPACT,DPAR) C C This is overdoing it for arrays with equal binwidth and pixel C filling factors of 1, (i.e. next run through loop XR could be C copied into the new XL), but we want to be fancy (general). C C Now search for the XI(j's) falling between XL and XR C Recall that we set up everything in such a way that we are C increasing in the value of XI and that the first XL will fall C after XI(LL) and the last XR before XI(LH). We will update C LL each time (we can do so because XO and XI were assumed to be C monotonic (or at least constant, but that doesn't make sense C at all). C C Fast loop here to find full pixels falling into XL,XR C A complete pix is defined GE left, LT right border C C DL is the distance between the lower border of pix XI(j) and C XL - if positive or 0 then we found the first applic. pixel. C C found the pix DO 60, J = LL,LH,LIP C this is the low pix DL = XL - XI(J) - DLIP*WI(J) IF (DL.LT.0.0D0) THEN ML = J C all in one pix DH = XI(J) + DLIP*WI(J) - XH C goto one-pix handling IF (DH.GE.0.0D0) THEN MH = ML GO TO 110 ENDIF GOTO 70 ENDIF 60 CONTINUE C CALL STETER(41,' ERROR - no lower bound of pixel') C C DR is the distance between the left border of pix XI(j) and C XR - if positive or 0 then we found the last applic. pixel. C Note that ML+1 is always .LE.LH because of boundary search C C dist. to right border 70 DO 80 J = ML + 1,LH,LIP MH = J DH = XI(J) + DLIP*WI(J) - XH C JMP to MH integration IF (DH.GT.0.0D0) THEN GO TO 90 ELSE C JMP to ML integration VAL = VAL + SNGL(YI(J)) C instead of : IF (DH.EQ.0.0D0) GO TO 100 IF (DABS(DH) .LT. 1.D-30) GOTO 100 ENDIF 80 CONTINUE C CALL STETER(41,' ERROR - no upper bound of pixel') C C handle the integrations (MH.NE.ML), in MH first C 90 PL = XI(MH) - DLIP*WI(MH) XINCR = (XH-PL)/XVLS CALL REBPIX(IOPT,NIN,XI,YI,WI,PL,XH,MINT,XINCR,MH,LIP,VALI) VAL = VAL + VALI/WI(MH) C C and look into ML in any case C 100 DL = XL - XI(ML) + DLIP*WI(ML) IF (DL.EQ.0.0D0) THEN VAL = VAL + SNGL(YI(ML)) ELSE PH = XI(ML) + DLIP*WI(ML) XINCR = (PH-XL)/XVLS CALL REBPIX(IOPT,NIN,XI,YI,WI,XL,PH,MINT,XINCR,ML,LIP, + VALI) VAL = VAL + VALI/WI(ML) END IF C GO TO 120 C C handle the single pix case (MH.EQ.ML) C C only one full pixel 110 DL = XL - XI(ML) + DLIP*WI(ML) IF (DL.EQ.0.0D0 .AND. DH.EQ.0.0D0) THEN VAL = SNGL(YI(MH)) ELSE XINCR = (XH-XL)/XVLS CALL REBPIX(IOPT,NIN,XI,YI,WI,XL,XH,MINT,XINCR,ML,LIP, + VALI) VAL = VALI/WI(ML) END IF C C update cuts and result array, CONSERVE FLUX IN PIXEL DOMAIN C 120 YO(I) = VAL IF (I.EQ.1) THEN RMIN = YO(1) RMAX = RMIN ELSE RMIN = AMIN1(RMIN,YO(I)) RMAX = AMAX1(RMAX,YO(I)) ENDIF C C update the loop start LL C LL = MH 130 CONTINUE C C Here is the end of REBMET C RETURN END C SUBROUTINE REBPIX(JOP,NP,DX,DY,WI,XA,XB,NPOI,DINCR,ISTA,IDIR,F) C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C Determine flux in fractions of pixel using various options JOP C C.INPUT C JOP I*4 option C NP I*4 number of array elements DX,DY C DX R*8 independent array C DY R*8 dependent array C WI R*8 bin (or pixel) widths C XA R*8 lower boundary of area of interest C XB R*8 upper boundary C NPOI I*4 number of points for integration/interpolation C to be used C DINCR R*8 interval between succesive points NPOI C ISTA I*4 start index or search loop in DX C IDIR I*4 Direction of DX array (X(1) < X(2) => -1) C C.OUTPUT C F R*4 intensity between XA and XB under DY = f(DX) C C ------------------------------------------------------------------- C INTEGER JOP,NP,NPOI,ISTA,IDIR REAL F DOUBLE PRECISION DX(NP),DY(NP),WI(NP),XA,XB DOUBLE PRECISION DINCR C GO TO (10,20,30),JOP CD WRITE (6,*) ' The option requested for interp./integ. is not' CD. 'possible' CD WRITE (6,*) ' Simple proportion of pixel will be used' GO TO 20 C C Spline interpolation and gaussian integration C 10 CALL REBISP(NP,DX,DY,XA,XB,NPOI,DINCR,ISTA,IDIR,F) RETURN C C Intensity is simply the corresponding proportion of the pixel's C intensity 20 F = DY(ISTA)* (XB-XA) RETURN C C linear interpolation and integration 30 CALL REBLNR(NP,DX,DY,XA,XB,ISTA,IDIR,F) RETURN C END SUBROUTINE REBISP(NP,DX,DY,XA,XB,NPOI,DINCR,IST,IDI,F) C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C F. Murtagh ST-ECF Version 1.0 15-AUG-1985 C M. Rosa ST-ECF 2.0 15-OCT-1985 C C Set up for REBIPL and integrate --> F C INTEGER I,NP,NPOI,IST,IDI,IFAIL,IS,ID,FLAG REAL F DOUBLE PRECISION DX(NP),DY(NP),XA,XB DOUBLE PRECISION VECI(10),VECD(10) DOUBLE PRECISION DINCR,DVAL,DER,XP C IFAIL = 0 IS = IST ID = IDI C XP = XA FLAG=1 CALL REBIPL(FLAG,XP,DVAL,DX,DY,NP,IS,ID) VECI(1) = XP VECD(1) = DVAL C DO 10 I = 2,NPOI XP = XP + DINCR FLAG=0 CALL REBIPL(FLAG,XP,DVAL,DX,DY,NP,IS,ID) VECI(I) = XP VECD(I) = DVAL 10 CONTINUE C DVAL = 0 DER = 0 CALL FINDIFF(VECI,VECD,NPOI,DVAL,DER,IFAIL) F = SNGL(DVAL) RETURN END SUBROUTINE REBIPL(FLAG,XP,P,X,F,N,IS,ID) C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C Interpolate a function P at a given value XP using C the table of values (X,F). C Spline interpolation scheme based on Hermite polynomials. C C Ref.: U.S. Airforce Surveys in Geophysics, no.272 C Publ of the Dominion Astrophys. Obs. 204, 1982 C C.AUTHOR J.D. PONZ (ESO) C C.MODIFIED M. ROSA (STECF) 02-OCT-1985 C C DOUBLE PRECISION calculation C Stepping through X is always increasing X and XP C All checking for validity range done outside C Loop parameters IS (start) and ID (direction) C supplied from outside C C Arguments : C XP Input argument C P Output interpolated value C X Input vector of independent values C F Input vector of dependent values C N Input number of points in (X,F) vector C IS Input index of respective pixel C ID Input direction of X array C INTEGER N,IS,ID,KA,KZ,IRUN,I,FLAG DOUBLE PRECISION F(N),X(N) DOUBLE PRECISION XP,P DOUBLE PRECISION LP1,LP2,L1,L2,FP1,FP2,XPI1,XPI,L1S,L2S SAVE KA,KZ,IRUN,FP1,FP2,LP1,LP2 C IF (FLAG.EQ.1) THEN IF (ID.GT.0) THEN KA = MAX0(IS-ID,1) KZ = N ELSE KA = MIN0(IS-ID,N) KZ = 1 END IF IRUN = 0 ENDIF C C DO 10 I = KA,KZ,ID IF (XP.LT.X(I)) GO TO 20 10 CONTINUE P = 10.**36 RETURN C 20 I = I - ID IF (I.EQ.KA-1 .AND. IRUN.NE.0) GO TO 30 KA = I + ID IRUN = 1 LP1 = 1./ (X(I)-X(I+1)) LP2 = 1./ (X(I+1)-X(I)) C IF (I.EQ.1) THEN FP1 = (F(2)-F(1))/ (X(2)-X(1)) ELSE FP1 = (F(I+1)-F(I-1))/ (X(I+1)-X(I-1)) END IF C IF (I.EQ.N-1) THEN FP2 = (F(N)-F(N-1))/ (X(N)-X(N-1)) ELSE FP2 = (F(I+2)-F(I))/ (X(I+2)-X(I)) END IF C 30 XPI1 = XP - X(I+1) XPI = XP - X(I) L1 = XPI1*LP1 L2 = XPI*LP2 L1S = L1*L1 L2S = L2*L2 P = F(I)* (1.D0-2.D0*LP1*XPI)*L1S + + F(I+1)* (1.D0-2.D0*LP2*XPI1)*L2S + FP2*XPI1*L2S + + FP1*XPI*L1S C RETURN END SUBROUTINE REBFNC(NF,D,NPA,NPAC,PA) C++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C.AUTHOR M.ROSA (STECF) VERSION 1.0 02-OCT-1985 C F.MURTAGH (STECF) 1.1 11-NOV-1985 C (more options, check on existence) C C Dummy routine to call the requested function f[NF] C Returns D[new] as f(D[old]) C C Functions: 1 LIN D --> a + b*D C 2 POL D --> a + b*D + c*DD + ... C 3 INV D --> a + b/D C 4 EXP D --> a + b*exp(D*c) C 5 DEX D --> a + b*10.**(D*c) C 6 LOG D --> a + b*LOG(D*c) C 7 DLG D --> a + b*LOG10(D*c) C 8 IPO D --> a + b/(D-q) + c/((D-q)**2) + ... C where q = PA(NPAC) and a,b,c... = PA(1,2,...(NPAC-1)) C 9 U01 User defined function REBU01 C INTEGER NF,NPA,NPAC DOUBLE PRECISION D,PA(NPA) C GO TO (10,20,30,40,50,60,70,80,90) NF C WRITE (*,*) ' Function not existent, sorry - look up manual ' D = 0. RETURN C 10 CALL REBLIN(D,PA,NPA,NPAC) RETURN 20 CALL REBPOL(D,PA,NPA,NPAC) RETURN 30 CALL REBINV(D,PA,NPA,NPAC) RETURN 40 CALL REBEXP(D,PA,NPA,NPAC) RETURN 50 CALL REBDEX(D,PA,NPA,NPAC) RETURN 60 CALL REBLOG(D,PA,NPA,NPAC) RETURN 70 CALL REBDLG(D,PA,NPA,NPAC) RETURN 80 CALL REBIPO(D,PA,NPA,NPAC) RETURN 90 CALL REBU01(D,PA,NPA,NPAC) RETURN C END SUBROUTINE REBLIN(DV,P,NP,NPC) INTEGER NP,NPC DOUBLE PRECISION DV,P(NP) DV = P(1) + P(2)*DV RETURN END C SUBROUTINE REBPOL(DV,P,NP,NPC) INTEGER NP,NPC,I DOUBLE PRECISION DV,P(NP),VD !MP 900117 VD = 0.0D0 DO 10 I = NPC,1,-1 VD = VD*DV + P(I) 10 CONTINUE DV = VD RETURN END C SUBROUTINE REBINV(DV,P,NP,NPC) INTEGER NP,NPC DOUBLE PRECISION DV,P(NP) DV = P(1) + P(2)/DV RETURN END C SUBROUTINE REBEXP(DV,P,NP,NPC) INTEGER NP,NPC DOUBLE PRECISION DV,P(NP) DV = P(1) + P(2)*EXP(DV*P(3)) RETURN END C SUBROUTINE REBDEX(DV,P,NP,NPC) INTEGER NP,NPC DOUBLE PRECISION DV,P(NP) DV = P(1) + P(2)*10.0D0** (DV*P(3)) RETURN END C SUBROUTINE REBLOG(DV,P,NP,NPC) INTEGER NP,NPC DOUBLE PRECISION DV,P(NP) DV = P(1) + P(2)*LOG(DV*P(3)) RETURN END C SUBROUTINE REBDLG(DV,P,NP,NPC) INTEGER NP,NPC DOUBLE PRECISION DV,P(NP) DV = P(1) + P(2)*LOG10(DV*P(3)) RETURN END C SUBROUTINE REBIPO(DV,P,NP,NPC) C C INVERSE POLYNOMIAL WITH OFFSET C Implement the PRISM formula for dispersion relations C C y = a + bz + cz**2 + dz**3 + .... P(npc-1)z**(npc-1) C C where z = x - p(npc) C INTEGER NP,NPC,I DOUBLE PRECISION P(NP),DV,VD VD = 0.0D0 DV = 1./ (DV-P(NPC)) IF (NPC.EQ.1) RETURN DO 10 I = NPC - 1,1,-1 VD = VD*DV + P(I) 10 CONTINUE DV = VD RETURN END C SUBROUTINE REBLNR(NP,DX,DY,XA,XB,ISTA,IDIR,F) C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C Subroutine linear C C F. Murtagh ST-ECF Version 1.0 17-Oct-1985 C C linear interpolation and subsequent integration. C C.INPUT C NP I*4 Number of array elements DX, DY C DX R*8 Independent array C DY R*8 Dependent array C XA R*8 Lower boundary of area of interest C XB R*8 Upper boundary C ISTA I*4 Start index for search loop in DX; C index number of pixel in which both C XA and XB are located C IDIR I*4 Direction of DX array (X(1) < X(2) => -1) C C.OUTPUT C F R*4 Intensity between XA and XB under DY = f(DX) C C------------------------------------------------------------------ INTEGER NP,ISTA,IDIR REAL F,YA,YB,XLIN DOUBLE PRECISION DX(NP),DY(NP),XA,XB DOUBLE PRECISION X,Y,X2,Y2 C C For convenience: X = DX(ISTA) Y = DY(ISTA) C C Four cases arise: XA or XB = X; otherwise, XA and XB are on the same C side of the pixel-centre, or they are on different sides. IF (XA.EQ.X) THEN IF (XB-X.GE.0) THEN Y2 = DY(ISTA+IDIR) X2 = DX(ISTA+IDIR) ELSE Y2 = DY(ISTA-IDIR) X2 = DX(ISTA-IDIR) END IF YB = XLIN(Y2,Y,X,XB,X2) F = ABS(X-XB)* (Y+YB)/2. ELSE IF (XB.EQ.X) THEN IF (XA-X.GE.0) THEN Y2 = DY(ISTA+IDIR) X2 = DX(ISTA+IDIR) ELSE Y2 = DY(ISTA-IDIR) X2 = DX(ISTA-IDIR) END IF YA = XLIN(Y2,Y,X,XA,X2) F = ABS(X-XA)* (Y+YA)/2. ELSE IF (((XA-X)/ (XB-X)).GT.0) THEN C (If XA and XB are both on the same side of the pixel-centre, C we get num. and denom. both of same sign, hence above is true: IF (XA-X.GT.0) THEN Y2 = DY(ISTA+IDIR) X2 = DX(ISTA+IDIR) ELSE Y2 = DY(ISTA-IDIR) X2 = DX(ISTA-IDIR) END IF YA = XLIN(Y2,Y,X,XA,X2) YB = XLIN(Y2,Y,X,XB,X2) F = (XB-XA)* (YA+YB)/2. ELSE IF (((XA-X)/ (XB-X)).LE.0) THEN IF (XA-X.GT.0) THEN YA = XLIN(DY(ISTA+IDIR),Y,X,XA,DX(ISTA+IDIR)) YB = XLIN(DY(ISTA-IDIR),Y,X,XB,DX(ISTA-IDIR)) ELSE YA = XLIN(DY(ISTA-IDIR),Y,X,XA,DX(ISTA-IDIR)) YB = XLIN(DY(ISTA+IDIR),Y,X,XB,DX(ISTA+IDIR)) END IF F = ABS(X-XA)* (Y+YA)/2. + ABS(X-XB)* (Y+YB)/2. END IF C RETURN END REAL FUNCTION XLIN(Y2,Y1,X1,X,X2) C C Function which returns linearly interpolated value corresponding to C X, where we have X1 <= X <= X2 or vice versa, and Y1 = f(X1) and C Y2 = f(X2). C DOUBLE PRECISION Y1,Y2,X,X1,X2 C XLIN = (Y2-Y1)* (X-X1)/ (X2-X1) + Y1 C END SUBROUTINE FINDIFF(X,Y,N,INT,ERR,IFAIL) C This subroutine integrates a function which is specified numerically C at four or more points, using third-order finite-difference formulae C according to a method due to Gill and Miller. (Computer Journal 15, C page 80-83,1972). C Author: Michele Peron, ESO-IPG. August 1990 INTEGER N,I,J,K,IFAIL DOUBLE PRECISION X(N),Y(N),INT,ERR,E,H1,H2,H3,H4,R1,R2,R3,R4 DOUBLE PRECISION D1,D2,D3,C,S,H3S IFAIL = 0 J = 3 K = N-1 ERR = 0 E = 0 INT=0 H1 = 0 H2 = 0 H3 = 0 H4 = 0 R1 = 0 R2 =0 R3 = 0 R4 = 0 D1 = 0 D2 = 0 D3 = 0 C = 0 S = 0 DO 10 I=J,K IF (I.EQ.J) THEN H2 = X(J-1)-X(J-2) D3 = (Y(J-1)-Y(J-2))/H2 H3 = X(J)-X(J-1) D1 = (Y(J)-Y(J-1))/H3 H1 = H2+H3 D2 = (D1-D3)/H1 H4 = X(J+1)-X(J) R1 = (Y(J+1)-Y(J))/H4 R2 = (R1-D1)/(H4+H3) H1 = H1+H4 R3 = (R2-D2)/H1 INT = H2*(Y(1)+H2*(D3/2.0-H2* 1 (D2/6.-(H2+2.*H3)*R3/12.))) S = -H2**3*(H2*(3.*H2+5.*H4)+10.*H3*H1)/60. ELSE H4 = X(I+1)-X(I) R1 = (Y(I+1)-Y(I))/H4 R4 = H4+H3 R2 = (R1-D1)/R4 R4 = R4+H2 R3 = (R2-D2)/R4 R4 = R4+H1 R4 = (R3-D3)/R4 ENDIF H3S = H3**2 INT = INT+H3*((Y(I)+Y(I-1))/2.-H3S* 1 (D2+R2+(H2-H4)*R3)/12.) C = H3*H3S*(2.*H3S+5.*(H3*(H4+H2)+ 2 2.*H4*H2))/120. E = E+(C+S)*R4 IF (I.EQ.J) THEN S = 2*C+S ELSE S = C ENDIF IF (I.EQ.K) THEN INT = INT+H4*(Y(N)-H4*(R1/2.+H4*(R2/6.+ 1 (2.*H3+H4)*R3/12.))) E = E-H4**3*R4*(H4*(3.*H4+5.*H2)+10. 2 *H3*(H2+H3+H4))/60. E = E+S*R4 ELSE H1 = H2 H2 = H3 H3 = H4 D1 = R1 D2 = R2 D3 = R3 ENDIF 10 CONTINUE ERR = E INT = INT+E END