|
|
SUBROUTINE BICUBSPL(X,Y,Z,NX,NY,NXDIM,NYDIM,AA,BB,CC)
Input: X,Y,Z,NX,NY,NXDIM,NYDIM
Output: AA,BB,CC
This subroutine computes the spline coefficients required by the subroutine
BICUBSPLX to compute a surface through bicubic spline interpolation.
REAL X(NX) -> X-values to be fitted
REAL Y(NY) -> Y-values to be fitted
REAL Z(NX,NY) -> Z-values to be fitted
INTEGER NX -> logical dimension of X and logical first dimension of Z
INTEGER NY -> logical dimension of Y and logical second dimension of Z
INTEGER NXDIM -> physical dimension of X and physical first dimension of Z
INTEGER NYDIM -> pyshical dimension of Y and physical second dimension of Z
REAL AA(NX,NY) -> spline coefficients of the one-dimensional cubic spline
fits to the rows of Z
REAL BB(NX,NY) -> spline coefficients of the one-dimensional cubic spline
fits to the rows of Z
REAL CC(NX,NY) -> spline coefficients of the one-dimensional cubic spline
fits to the rows of Z
SUBROUTINE BICUBSPLX(X,Y,Z,NX,NY,NXDIM,NYDIM,AA,BB,CC,X0,Y0,Z0)
Input: X,Y,Z,NX,NY,NXDIM,NYDIM,AA,BB,CC,X0,Y0
Output: Z0
This subroutine computes the bicubic spline Z0 at X0,Y0, using the spline
coefficients computed with BICUBSPL.
REAL X(NX) -> X-values fitted with BICUBSPL
REAL Y(NY) -> Y-values fitted with BICUBSPL
REAL Z(NX,NY) -> Z-values fitted with BICUBSPL
INTEGER NX -> logical dimension of X and logical first dimension of Z
INTEGER NY -> logical dimension of Y and logical second dimension of Z
INTEGER NXDIM -> physical dimension of X and physical first dimension of Z
INTEGER NYDIM -> pyshical dimension of Y and physical second dimension of Z
REAL AA(NX,NY) -> coefficients of the one-dimensional cubic spline fits
to the rows of Z computed with BICUBSPL
REAL BB(NX,NY) -> coefficients of the one-dimensional cubic spline fits
to the rows of Z computed with BICUBSPL
REAL CC(NX,NY) -> coefficients of the one-dimensional cubic spline fits
to the rows of Z computed with BICUBSPL
REAL X0 -> X-value where the bicubic spline will be evaluated
REAL Y0 -> Y-value where the bicubic spline will be evaluated
REAL Z0 -> bicubic spline value at X0,Y0
SUBROUTINE BINSEARCH(X,N,X0,N0)
Input: X,N,X0,N0 Output: N0 Given the array X(N), and the test value X0, this subroutine returns an integer N0, such that X0 is between X(N0) and X(N0+1). As input N0 is employed to start the searching. If X0.LT.X(1) then N0=0 on output, whereas if X0.GT.X(N) then N0=N. If X0.EQ.X(K), N0=K on output. REAL X(N) -> ordered input array (not necesarilly equally-spaced) INTEGER N -> no. of points in input array REAL X0 -> argument to be searched for INTEGER N0 -> location of X0 in the input array
SUBROUTINE BROADEN(S1,S2,NCHAN,STWV,DISP,SIGMA,LERR)
Input: S1,NCHAN,STWV,DISP,SIGMA,LERR
Output: S2
Broadens a single spectrum by convolving with a gaussian (variable along the
wavelength scale).
REAL S1(NCMAX) -> input spectrum
REAL S2(NCMAX) -> output spectrum
INTEGER NCHAN -> no. of channels
REAL STWV -> central wavelength of the first channel
REAL DISP -> dispersion (Angstrom/pixel)
REAL SIGMA(NCMAX) -> sigma value of the gaussian to be applied in each
channel
LOGICAL LERR -> if .TRUE. S1 corresponds to an error spectrum
SUBROUTINE CAUCHYFIT(X0,SIGMA,AMP,EEX0,EESIGMA,EEAMP,YRMSTOL)
Input: YRMSTOL Input (COMMON): NP,XF,YF Output: X0,SIGMA,AMP,EEX0,EESIGMA,EEAMP Fit numerically a Cauchy function (using DOWNHILL): Y=AMP/[SIGMA^2+(X-X0)^2] REAL X0 -> center of the fitted Cauchy function REAL SIGMA -> sigma value of the fitted Cauchy function REAL AMP -> maximum of the fitted Cauchy function REAL EEX0,EESIGMA,EEAMP -> errors in X0,SIGMA,AMP (rms from DOWNHILL) REAL YRMSTOL -> stopping criterion for DOWNHILL
SUBROUTINE CFFTD(N,XR,XI,IMODE)
Input N,XR,XI,IMODE Output XR,YR If IMODE=1, this subroutine computes the FFT of the input complex vector XR + i XI, where XR and XI are both real variables. As output XR and XI are the real and imaginary part of the transform. If IMODE=-1 the subroutine evaluates the inverse FFT. See E. O. Brigham, The Fast Fourier Transform, pag.160. INTEGER N -> number of points (must be a power of 2) DOUBLE PRECISION XR(N) -> real part of the input data DOUBLE PRECISION XI(N) -> imaginary part of the input data INTEGER IMODE -> +1 (direct FFT) or -1 (inverse FFT)
SUBROUTINE CFFT(N,XR,XI,IMODE)
Input N,XR,XI,IMODE Output XR,YR If IMODE=1, this subroutine computes the FFT of the input complex vector XR + i XI, where XR and XI are both real variables. As output XR and XI are the real and imaginary part of the transform. If IMODE=-1 the subroutine evaluates the inverse FFT. See E. O. Brigham, The Fast Fourier Transform, pag.160. INTEGER N -> number of points (must be a power of 2) REAL XR(N) -> real part of the input data REAL XI(N) -> imaginary part of the input data INTEGER IMODE -> +1 (direct FFT) or -1 (inverse FFT)
SUBROUTINE CHEQUEA_FILEINDEX
This subroutine verifies that there is an appropiate file containing the index definitions.
SUBROUTINE CHREBIN(CHANSHIFT,NCHAN,S,SS)
Input: CHANSHIFT,NCHAN,S Output: SS This subroutine applies a constant channel shift to a spectrum. REAL CHANSHIFT -> channel shift to be applied INTEGER NCHAN -> number of channels REAL S(NCHAN) -> initial spectrum REAL SS(NCHAN) -> shifted spectrum
DOUBLE PRECISION FUNCTION COMBPF(N,K)
Input: N,K Output: COMBPF (function) Calculate the binomial coefficient N over K INTEGER N INTEGER K
SUBROUTINE CUBSPL(X,Y,N,IMODE,S,A,B,C)
Input: X,Y,N,IMODE,S
Output: A,B,C,S
This subroutine computes the coefficients of a cubic spline. See C.F. Gerald
and P. O. Wheatley, in Applied Numerical Analysis, 4th edition, pag. 207.
The subroutine returns the spline coefficients, where the spline defined
in the interval between X(I),Y(I) and X(I+1),Y(I+1) is given by:
Y = A(I)*(X-X(I))**3 + B(I)*(X-X(I))**2 + C(I)*(X-X(I)) + D(I)
REAL X(N) -> X-values to be fitted
REAL Y(N) -> Y-values to be fitted
INTEGER N -> number of data points
INTEGER IMODE -> End conditions mode: if S(I) represent the second derivative
at the point X(I),Y(I), the following four possibilites
are available:
1) IMODE=1: S(0)=0, S(N)=0. This is called natural cubic
spline. It is equivalent to assuming that the end cubics
aproach linearity at their extremities.
2) IMODE=2: S(0)=S(1), S(N)=S(N-1). This is equivalent to
assuming that the cubics approach parabolas at their
extremities.
3) IMODE=3: S(0) is a linear extrapolation from S(1) and
S(2), and S(N) is a linear extrapolation from S(N-2)
and S(N-1).
4) IMODE=4: Force the slopes at each end to assume certain
values.
REAL S(N) -> if IMODE=4, in input S(1) and S(N) contain the first derivatives
at X(1) and X(N). In output, this matrix contains the second
derivatives
REAL A(N) -> spline coefficients
REAL B(N) -> spline coefficients
REAL C(N) -> spline coefficients
SUBROUTINE CUBSPLX(X,Y,A,B,C,N,I0,X0,Y0)
Input: X,Y,A,B,C,N,I0,X0
Output: Y0
The subroutine returns the cubic spline evaluated at X0, using the
coefficients determined in a previous call to CUBSPL. The spline defined in
the interval between X(I),Y(I) and X(I+1),Y(I+1) is given by:
Y = A(I)*(X-X(I))**3 + B(I)*(X-X(I))**2 + C(I)*(X-X(I)) + D(I)
If X0.LT.X(1), I=1 is employed (first computed spline)
If X0.GT.X(N), I=N-1 is employed (last computed spline)
REAL X(N) -> X-values fitted with CUBSPL
REAL Y(N) -> Y-values fitted with CUBSPL
REAL A(N) -> spline coefficients
REAL B(N) -> spline coefficients
REAL C(N) -> spline coefficients
INTEGER N -> number of data points
INTEGER I0 -> initial location to start the search of the place of X0 in
the X array
REAL X0 -> X-value where the spline function will be evaluated
REAL Y0 -> spline value at X0
SUBROUTINE DOWNHILL(N,X0,DX0,YFUNK,A,B,G,YRMSTOL,XF,DXF,NEVAL)
Input N,X0,DX0,YFUNK,A,B,G,YRMSTOL
Output XF,NEVAL
Minimization of the function YFUNK of N variables using the downhill
simplex method, as explained by Nelder and Mead (1965, Computer Journal, 7,
pags. 308-313). The routine returns when the stopping criterion is reached
(the r.m.s. of the YFUNK values computed with all the vertices of the simplex
is .LT. YRMSTOL), or when the number of function evaluations
is too large (NEVAL.GT.NEVALMAX).
INTEGER N -> number of variables
REAL X0(N) -> starting point (initial solution)
REAL DX0(N) -> N characteristic length scales, employed to derive N
aditional starting points which, together with X0, form
the (N+1) vertices of the simplex
REAL YFUNK -> function to be minimized
REAL A -> reflection coefficient (ALPHA); tipically, A=1.0
REAL B -> contraction coefficient (BETA); tipically, B=0.5
REAL G -> expansion coefficient (GAMMA); tipically, G=2.0
REAL YRMSTOL -> stopping criterion
REAL XF(N) -> final solution
REAL DXF(N) -> rms of XF evaluated from the final different points of the
simplex
INTEGER NEVAL -> number of evaluations of YFUNK employed by DOWNHIL to reach
the solution; the routine returns NEVAL=-1 if something goes
wrong
DOUBLE PRECISION FUNCTION FACTORIALPF(N)
Input: N Output: FACTORIALPF (function) Calculate N factorial INTEGER N
SUBROUTINE FFT2POWER(N0,N)
Input: N0 Output: N Given an integer N0, this subroutine asks for a power of 2, such as N=2**K, being K integer and N.GE.N0. This value of N is employed by other subroutines to perform zero padding prior computing FFT.
SUBROUTINE FFTCORREL(N,DATA1,DATA2,XCORR,FCORR)
Input N,DATA1,DATA2
Output XCORR,FCORR
Compute the correlation function FCORR of two real data vectors DATA1 and
DATA2, using FFT. We assume that both data sets have been properly filtered
and dimensioned. We use the discrete correlation theorem, which says that
the discrete correlation of two real functions f1 and f2 is one member of
the discrete Fourier transform pair:
Corr(f1,f2) <==> F1 F2*
where F1 and F2 are the discrete Fourier transforms of f1 and f2,
respectively, and asterisk denotes complex conjugation.
Note that DATA1 and DATA2 are not modified by this routine.
INTEGER N -> number of points (must be a power of 2)
REAL DATA1(N) -> first data set to be correlated
REAL DATA2(N) -> second data set to be correlated
REAL XCORR(N) -> abcissas of the output correlation function
REAL FCORR(N) -> output correlation function
SUBROUTINE FFTCORRZOOM(X,Y,N,NL,NAME1,NAME2,NFIT,LPLOT,X0,Y0)
Input: X,Y,N,NL,NAME1,NAME2,NFIT,LPLOT
Output: X0,Y0
This subroutine plots the correlation function given by X(N),Y(N), zooms in
around the maximum and fit a second-order polynomial in order to obtain the
exact location of the peak. The fitted value is returned through X0,Y0.
REAL X(N) -> X-coordinates of the correlation function
REAL Y(N) -> Y-coordinates of the correlation function
INTEGER N -> dimension of X and N
INTEGER NL -> number of pixels affected by zero padding
CHARACTER*(*) NAME1 -> description of the first data set employed to compute
cross correlation
CHARACTER*(*) NAME2 -> description of the second data set employed to compute
cross correlation
INTEGER NFIT -> no. of points around peak to fit maximum (if NFIT=0, the
routine asks for this number; otherwise the routine returns
without prompting)
LOGICAL LPLOT -> if .TRUE., plot zoomed region
REAL X0 -> X-offset of the peak of the correlation function
REAL Y0 -> peak value of the correlation function
SUBROUTINE FFTCOSBELL(N,COSBELL,FL)
Input: N, FL Output: COSBELL Compute a cosine bell expanding N pixels, being FL the fraction of pixels (at the beginning and at the end of the cosine bell) employed to perform the transition from zero to one. See Brault & White, A&A, 13, 169. INTEGER N -> dimension of COSBELL (pixels) REAL COSBELL(N) -> cosine bell REAL FL -> fraction of pixels at the borders of the cosine bell
SUBROUTINE FFTKFILTER(N,KFILTER,K1,K2,K3,K4)
Input: N,K1,K2,K3,K4 Output: KFILTER Given K1, K2, K3 and K4, this subroutine computes the filter KFILTER, such as: KFILTER(I)=0., I=1,...,K1 KFILTER(I) is a straight line from 0. to 1., I=K1,...,K2 KFILTER(I)=1., I=K2,...,K3 KFILTER(I) is a straight line from 1. to 0., I=K3,...,K4 KFILTER(I)=0., I=K4,...,N/2+1 KFILTER(I)=0., I=N,...,N-K1+2 KFILTER(I) is a straight line from 0. to 1., I=N-K1+2,...,N-K2+2 KFILTER(I)=1., I=N-K2+2,...,N-K3+2 KFILTER(I) is a straight line from 1. to 0., I=N-K3+2,...,N-K4+2 KFILTER(I)=0., I=N-K4+2,...,N/2+1 N must be a power of 2, and 1.le.K1.le.K2.le.K3.le.K4.le.(N/2-1). This filter is employed by other routines to perform a frequency filter in the frequency domain of the FFT.
SUBROUTINE FFTPREP(N,SP,NC1,NC2,LNORM,COSBELL,LFILT,KFILTER,LPLOT,CNAME)
Input: N,SP,NC1,NC2,LNORM,COSBELL,LFILT,KFILTER,LPLOT,CNAME Output: SP This subroutine prepares spectrum SP for cross correlation. For this purpose the following steps are performed: - normalization of the spectrum, using only the data in the range [SP(NC1),SP(NC2)] - extraction of the useful spectral region: [SP(NC1),SP(NC2)] is converted into [SP(1),SP(NC2-NC1+1)] - multiplication of the spectrum by a cosine bell defined from 1 to NC2-NC1+1 - zero padding from NC2-NC1+2 to N, where N=2**K, being K an integer - computation of the Fast Fourier Transform of SP(1),...,SP(N) and obtention of the power spectrum - filtering of the power spectrum by applying the user defined filter KFILTER - computation of the inverse FFT of the filtered power spectrum and obtention of the filtered spectrum SP(1),...,SP(N) - extraction of the initial spectral region: SP(1),...,SP(NC2-NC1+1) - multiplication by the cosine bell defined from 1 to NC2-NC1+1 (again!) - zero padding from NC2-NC1+2 to N (again!) INTEGER N -> number of points in output (must be a power of 2) REAL SP(N) -> spectrum to be prepared for cross correlation INTEGER NC1 -> index which indicates the first SP value to be employed INTEGER NC2 -> index which indicates the last SP value to be employed LOGICAL LNORM -> if .TRUE. SP is normalized prior any other manipulation REAL COSBELL(N) -> cosine bell (defined from point number 1 to NC2-NC1+1) LOGICAL LFILT -> if .TRUE. SP is filtered applying the filter KFILTER REAL KFILTER(N) -> filter (defined from point number 1 to N) LOGICAL LPLOT -> if .TRUE. plot intermediate plots CHARACTER*(*) CNAME -> string description for plots
REAL FUNCTION FINTGAUSS(X1,X2,N,X0,FACTOR1) -> integral of Gaussian REAL FUNCTION FINTGAUSSE(X1,X2,N,X0,FACTOR1) -> integral of Gaussian^2
Input: X1,X2,N,X0,FACTOR1 Output: FINTGAUSS (function) or FINTGAUSSE (function) Calculate the integral of a Gaussian of width SIGMA, center at wavelength X0, using Simpson's rule, between wavelengths X1 and X2, with N intervals (N must be even). REAL X1,X2 -> integral limits (wavelengths) INTEGER N -> number of intervals to estimate integral REAL X0 -> Gaussian center (wavelength) REAL FACTOR1 -> = (c^2/(2*SIGMA^2))/(X0^2)
REAL FUNCTION FMEAN0(N,X,SIGMA)
Input: N,X Output: FMEAN0 (function),SIGMA Calculate the mean value of X(N) and its r.m.s. INTEGER N -> no. of elements REAL X(N) -> input matrix REAL SIGMA -> r.m.s. around the mean value
REAL FUNCTION FMEAN1(N,X)
Input: N,X Output: FMEAN1 (function) Calculate the mean value of X(N) INTEGER N -> no. of elements REAL X(N) -> input matrix
REAL FUNCTION FMEAN2(N,X,TIMES)
Input: N,X,TIMES Output: FMEAN2 (function) Calculate the mean value of X(N) rejecting points at TIMES sigma. The function can recover points rejected in previous iterations. INTEGER N -> no. of elements REAL X(N) -> input matrix REAL TIMES -> times sigma to reject points before calculating the mean
REAL FUNCTION FMEDIAN1(N,X)
Input: N,X Output: FMEDIAN (function), X (sorted) Calculate the median value of X(N). It is important to note that this subroutine rearranges the matrix X which is returned sorted. INTEGER N -> no. of elements REAL X(N) -> input matrix
REAL FUNCTION FPERCENT(N,X,PERCENTILE)
Input: N,X,PERCENTILE Output: FPERCENT (function) Calculate a fixed percentile of X(N) INTEGER N -> no. of elements REAL X(N) -> input matrix REAL PERCENTILE -> percentile to be computed
REAL FUNCTION FPOLY(NDEG,COEFF,X)
Input: NDEG,COEFF,X Output: FPOLY (function) Evaluate the polynomial of degree NDEG and coefficients COEFF at X. INTEGER NDEG -> polynomial degree REAL COEFF(NDEG+1) -> polynomial coefficients REAL X -> abscissa at which the polynomial is going to be evaluated
SUBROUTINE GAUSCFIT(X0,SIGMA,AMP,Y0,EEX0,EESIGMA,EEAMP,EEY0,YRMSTOL)
Input: YRMSTOL Input (COMMON): NP,XF,YF Output: X0,SIGMA,AMP,Y0,EEX0,EESIGMA,EEAM,EEY0 Fit numerically a gaussian + constant (using DOWNHILL): Y=Y0+AMP*EXP[-((X-X0)^2/(2*SIGMA^2))] REAL X0 -> center of the fitted gaussian REAL SIGMA -> sigma value of the fitted gaussian REAL AMP -> maximum of the fitted gaussian REAL Y0 -> constant REAL EEX0,EESIGMA,EEAMP,EEY0 -> errors in X0,SIGMA,AMP,Y0 (rms from DOWNHILL) REAL YRMSTOL -> stopping criterion for DOWNHILL
SUBROUTINE GAUSCFIT_MOVEL(X0,SIGMA,AMP,Y0,EEX0,EESIGMA,EEAMP,EEY0,YRMSTOL)
Input: YRMSTOL Input (COMMON): NP,XF,YF Output: X0,SIGMA,AMP,Y0,EEX0,EESIGMA,EEAM,EEY0 Fit numerically a gaussian + constant (using DOWNHILL): Y=Y0+AMP*EXP[-((X-X0)^2/(2*SIGMA^2))] REAL X0 -> center of the fitted gaussian REAL SIGMA -> sigma value of the fitted gaussian REAL AMP -> maximum of the fitted gaussian REAL Y0 -> constant REAL EEX0,EESIGMA,EEAMP,EEY0 -> errors in X0,SIGMA,AMP,Y0 (rms from DOWNHILL) REAL YRMSTOL -> stopping criterion for DOWNHILL This subroutine is identical to GAUSCFIT, but the physical dimensions of the data to be fitted is defined to be NMAXFFT. This change is required for the program movel.
SUBROUTINE GAUSS2AFIT(NPFIT,XFIT,YFIT,EYFIT,DELTAX,X0,SIGMA,AMP,
EX0,ESIGMA,EAMP,EEX0,EESIGMA,EEAMP,
YRMSTOL,NSIMUL)
Input: NPFIT,XFIT,YFIT,EYFIT,DELTAX,YRMSTOL,NSIMUL Output: X0,SIGMA,AMP,EX0,ESIGMA,EAMP,EEX0,EESIGMA,EEAMP Numerical fit of 2 gaussians with the same width and area (using DOWNHILL), with a fixed separation given by DELTAX. Y=AMP*EXP[-((X-X0)^2/(2*SIGMA^2))]+AMP*EXP[-((X-X0-DELTAX)^2/(2*SIGMA^2))] INTEGER NPFIT -> number of points to be fitted REAL XFIT,YFIT,EYFIT -> x, y and error REAL DELTAX -> separation between gaussians REAL X0 -> center of the fitted gaussian REAL SIGMA -> sigma value of the fitted gaussian REAL AMP -> maximum of the fitted gaussian REAL EX0,ESIGMA,EAMP -> errors in X0,SIGMA,AMP (due to EYFIT --simulations--) REAL EEX0,EESIGMA,EEAMP -> errors in X0,SIGMA,AMP (rms from DOWNHILL) REAL YRMSTOL -> stopping criterion for DOWNHILL INTEGER NSIMUL -> number of simulations to compute errors
SUBROUTINE GAUSS2BFIT(NPFIT,XFIT,YFIT,EYFIT,DELTAX,X0,SIGMA,AMP1,AMP2,
EX0,ESIGMA,EAMP1,EAMP2,EEX0,EESIGMA,EEAMP1,EEAMP2,
YRMSTOL,NSIMUL)
Input: NPFIT,XFIT,YFIT,EYFIT,DELTAX,YRMSTOL,NSIMUL Output: X0,SIGMA,AMP1,AMP2,EX0,ESIGMA,EAMP1,EAMP2,EEX0,EESIGMA,EEAMP1,EEAMP2 Numerical fit of 2 gaussians with the same width and different area (using DOWNHILL), with a fixed separation given by DELTAX. Y=AMP1*EXP[-((X-X0)^2/(2*SIGMA^2))]+AMP2*EXP[-((X-X0-DELTAX)^2/(2*SIGMA^2))] INTEGER NPFIT -> number of points to be fitted REAL XFIT,YFIT,EYFIT -> x, y and error REAL DELTAX -> separation between gaussians REAL X0 -> center of the fitted gaussian REAL SIGMA -> sigma value of the fitted gaussian REAL AMP1 -> maximum of the fitted gaussian #1 REAL AMP2 -> maximum of the fitted gaussian #2 REAL EX0,ESIGMA,EAMPn-> errors in X0,SIGMA,AMPn (due to EYFIT --simulations--) REAL EEX0,EESIGMA,EEAMPn-> errors in X0,SIGMA,AMPn (rms from DOWNHILL) REAL YRMSTOL -> stopping criterion for DOWNHILL INTEGER NSIMUL -> number of simulations to compute errors
SUBROUTINE GAUSSFITAMP(NPFIT,XFIT,YFIT,EYFIT,X0,SIGMA,EX0,ESIGMA,MINSIGMA,
AMP,EAMP,EEAMP,NSIMUL)
Input: NPFIT,XFIT,YFIT,EYFIT,X0,SIGMA,EX0,ESIGMA,MINSIGMA,NSIMUL
Output: AMP,EAMP,EEAMP
Fit of AMP of a gaussian with X0 and SIGMA fixed:
Y=AMP*EXP[-((X-X0)^2/(2*SIGMA^2))]
INTEGER NPFIT -> number of points to be fitted
REAL XFIT,YFIT,EYFIT -> x, y and error
REAL X0, EX0 -> center of the gaussian and its error
REAL SIGMA, ESIGMA -> sigma value of the gaussian and its error
REAL MINSIGMA -> minimum SIGMA allowed in simulations (typically MINSIGMA
must be the spectral resolution)
REAL AMP -> maximum of the fitted gaussian
REAL EAMP -> error in SIGMA (due to EYFIT)
REAL EEAMP -> error in SIGMA (due to EX0 and ESIGMA ---simulations---)
INTEGER NSIMUL -> number of simulations to compute EEAMP
SUBROUTINE GAUSSFIT(NPFIT,XFIT,YFIT,EYFIT,X0,SIGMA,AMP,
EX0,ESIGMA,EAMP,EEX0,EESIGMA,EEAMP,YRMSTOL,NSIMUL)
Input: NPFIT,XFIT,YFIT,EYFIT,YRMSTOL,NSIMUL Output: X0,SIGMA,AMP,EX0,ESIGMA,EAMP,EEX0,EESIGMA,EEAMP Numerical fit of a gaussian (using DOWNHILL): Y=AMP*EXP[-((X-X0)^2/(2*SIGMA^2))] INTEGER NPFIT -> number of points to be fitted REAL XFIT,YFIT,EYFIT -> x, y and error REAL X0 -> center of the fitted gaussian REAL SIGMA -> sigma value of the fitted gaussian REAL AMP -> maximum of the fitted gaussian REAL EX0,ESIGMA,EAMP -> errors in X0,SIGMA,AMP (due to EYFIT --simulations--) REAL EEX0,EESIGMA,EEAMP -> errors in X0,SIGMA,AMP (rms from DOWNHILL) REAL YRMSTOL -> stopping criterion for DOWNHILL INTEGER NSIMUL -> number of simulations to compute errors
REAL FUNCTION INTEGTAB(N,X,Y,X1,X2,IFLAG1,IFLAG2)
Input: N,X,Y,X1,X2
Output: INTEGTAB(function), IFLAG1,IFLAG2
Performs the integration of a function given in a tabular form, between the
limits X1 and X2. Note that the X matrix must be sorted in ascending order.
INTEGER N -> input number of data in X and Y
REAL X(N) -> data matrix
REAL Y(N) -> data matrix
REAL X1 -> first limit of the integral
REAL X2 -> second limit of the integral
INTEGER IFLAG1 -> = 0 : interpolation
= -1 : extrapolation towards lower X values
= +1 : extrapolation towards higher X values
= +9 : error (division by zero)
INTEGER IFLAG2 -> = 0 : interpolation
= -1 : extrapolation towards lower X values
= +1 : extrapolation towards higher X values
= +9 : error (division by zero)
SUBROUTINE LAGRANGE(N,X,Y,X0,Y0)
Input: N,X,Y,X0 Output: Y0 Calculate the Lagrangian polynomial and evaluate such polynomial at X=X0. We do not assume uniform spacing between the x-values, nor do we need the x-values arranged in a particular order. However, the x-values must all be distinct. We follow the algorithm described by B.P. Demidovich and I.A. Maron in Calculo Numerico Fundamental, Paraninfo 1988, pag. 593. INTEGER N -> no. of elements REAL X(N) -> input x-matrix REAL Y(N) -> input y-matrix REAL X0 -> x-value where the Lagrangian polynomial is evaluated REAL Y0 -> polynomial value at X=X0
SUBROUTINE FITL(X,Y,NX,sig,iw,xmn,xmx,A,B,SIGA,SIGB,STD)
-----------------------------------------------------------------------C
Linear Fit Routines J. Jesus Gonzalez G.
-----------------------------------------------------------------------C
-----------------------------------------------------------------------C
Fits the line y = b*x + a
THIS ROUTINE ITERATES ELIMINATING HIGHLY DEVIANT POINTS
INPUT: X, Y - Data arrays
SIG - Y-error of points (<=0 if a point is to be rejected)
IW - =0 if unweighted fit, weighted fit otherwise.
xmn,xmx - Limits of fit.
OUTPUT:
B, A - Slope and zero-ordinate.
SIGB, SIGA - Stimated errors on B and A.
STD - Unbiased Std-Deviation from the fit.
Remeber how to compute errors of predicted values:
VAR(y(x)) = VAR(a) + VAR(b)*(x^2 - 2*xm*x), since the
ab-error covariance is COV(ab)=-xm*VAR(b)
-----------------------------------------------------------------------C
REAL FUNCTION LININTERP(N,X,Y,X0,IFLAG,N1,N2)
Input: N,X,Y,X0
Output: LININTERP(function), IFLAG
Performs a linear interpolation in the table X(N),Y(N) at x=X0. Note that the
X matrix must be sorted in ascending order, although the
INTEGER N -> input number of data in X and Y
REAL X(N) -> data matrix
REAL Y(N) -> data matrix
REAL X0 -> x-point at which the linear interpolation is evaluated
INTEGER IFLAG -> = 0 : interpolation
= -1 : extrapolation towards lower X values
= +1 : extrapolation towards higher X values
= +2 : X0=X(N), which could produce some "border" effects
= +9 : error (division by zero)
INTEGER N1 -> first data entry towards the left of X0
INTEGER N2 -> first data entry towards the right of X0
SUBROUTINE LUDCMP(A,N,NDIM,ORDER,SCALEROW,IOK,IPAR)
Input: A,N,NDIM
Output: A,ORDER,SCALEROW,IOK,IPAR
This subroutine computes the L and U triangular matrices equivalent to the
A matrix, such that LU = A. These matrices are returned in the space of A,
in compact form. See C.F. Gerald and P. O. Wheatley, in Applied Numerical
Analysis, 4th edition, pag. 106.
REAL A(NDIM,NDIM) -> matrix of coefficients
INTEGER N -> logical dimension of A
INTEGER NDIM -> physical dimension of A in the calling program
INTEGER ORDER(N) -> vector holding row order after pivoting
REAL SCALEROW(N) -> vector holding scaling factors applied to each row
INTEGER IOK -> returns 0 if everything works properly, +(the row number)
if all elements in a row are zero, or -(the row number) if
the pivot value is zero.
INTEGER IPAR -> returns as +1 or -1 depending on whether the number of row
interchanges was even or odd, respectively
SUBROUTINE LUSOLV(A,N,NDIM,ORDER,SCALEROW,B,X)
Input: A,N,NDIM,ORDER,SCALEROW,B Output: X This subroutine solves the set of N linear equations A X = B, where the A matrix corresponds to the LU decomposition of the initial coefficient matrix. See C.F. Gerald and P. O. Wheatley, in Applied Numerical Analysis, 4th edition, pag. 110. The matrix A remains unchanged (also ORDER and SCALEROW), so subsequent calls to this subroutine, variying the B matrix, can be performed. REAL A(NDIM,NDIM) -> matrix of coefficients (LU in compact scheme) INTEGER N -> logical dimension of A INTEGER NDIM -> physical dimension of A in the calling program INTEGER ORDER(N) -> vector holding row order after pivoting in LUDCMP REAL SCALEROW(N) -> vector holding scaling factors applied to each row REAL B(N) -> right-hand side vector B REAL X(N) -> solution vector X
INTEGER FUNCTION MIDEIND(NS1,NS2,ITI,WV,FWV,CERR,RCVEL1,NCRES,FINDEX,EINDEX,
EJJGG,ESIMU,SN,FFPLAW,LONLY_SN)
Input: NS1,NS2,ITI,WV,CERR,RCVEL1,NCRES,FFPLAW,LONLY_SN
Output: FINDEX,EINDEX,EJJGG,ESIMU,SN
Input/Output: (see COMMON blocks)
Return MIDEIND=0 if atomic/molecular/D4000/B4000/generic index (and error)
has been properly measured. MIDEIND=1 otherwise.
INTEGER NS1,NS2 -> scans to be coadded
INTEGER ITI -> index type:
ITI= 1: molecular index
ITI= 2: atomic index
ITI= 3: D4000 (Bruzual 1983)
ITI= 4: B4000 (own defintion)
ITI = 5: like B4000 but computing flux per angstrom
ITI=????: generic index:
ITI= C x 100 + L, where
C: number of continuum regions
L: number of absorption regions
Since Cmin=1, Cmax=99, Lmin=1, Lmax=99
==> ITImin=101, ITImax=9999
REAL WV(NWVMAX) -> wavelength limits
REAL FWV(NWVMAX/4) -> constant factors to be used when computing
the index as multiplicative coefficients for
the absorption signal.
CHARACTER*1 CERR -> 'y' if error spectrum is available, 'n' otherwise
REAL RCVEL1 -> 1 + v/c
INTEGER NCRES -> No. of response curve to be employed (1=averaged)
REAL FINDEX -> measured index
REAL EINDEX -> measured index error using own formulae
REAL EJJGG -> measured index error using JJGGs formulae
REAL ESIMU -> measured index error using numerical simulations
REAL SN -> averaged (S/N ratio)/Angstrom in the index region
REAL FFPLAW -> =0,...,10 indicates the fraction of light (in the
selected photometric band) to be used with a power law
of the form: F(lambda) = k * (lambda)**(alpha-2.0),
where alpha is defined through the variable ALPHAPLAW
(global variable ---COMMON---). If FFPLAW = -1 no
power law is employed.
LOGICAL LONLY_SN -> if .TRUE., the subroutine only computes the
mean S/N ratio and returns
SUBROUTINE ORDENA1F1I(N,A,B)
Input: N,A,B Output: A,B Sorts the real array A(N) into ascending numerical order. The additional array B(N) is simultaneously changed in parallel with the array A(N). Note that the two input arrays are returned rearranged. We follow the Heapsort method, as described by D. Knuth, The Art of Computer Programming (pag.146, 5.2.3.). INTEGER N -> input number of data in A REAL A(N) -> data matrix to be sorted INTEGER B(N) -> data matrix to be sorted in parallel with matrix A
SUBROUTINE ORDENA1F2I(N,A,B,C)
Input: N,A,B,C Output: A,B,C Sorts the real array A(N) into ascending numerical order. The additional arrays B(N) and C(N) are simultaneously changed in parallel with the array A(N). Note that the three input arrays are returned rearranged. We follow the Heapsort method, as described by D. Knuth, The Art of Computer Programming (pag.146, 5.2.3.). INTEGER N -> input number of data in A REAL A(N) -> data matrix to be sorted INTEGER B(N) -> data matrix to be sorted in parallel with matrix A INTEGER C(N) -> data matrix to be sorted in parallel with matrix A
SUBROUTINE ORDENA1F(N,A)
Input: N,A Output: A Sorts the real array A(N) into ascending numerical order. Note that the input array is returned rearranged. We follow the Heapsort method, as described by D. Knuth, The Art of Computer Programming (pag.146, 5.2.3.). INTEGER N -> input number of data in A REAL A(N) -> data matrix to be sorted
SUBROUTINE ORDENA1I(N,A)
Input: N,A Output: A Sorts the integer array A(N) into ascending numerical order. Note that the input array is returned rearranged. We follow the Heapsort method, as described by C D. Knuth, The Art of Computer Programming (pag.146, 5.2.3.). INTEGER N -> input number of data in A INTEGER A(N) -> data matrix to be sorted
SUBROUTINE ORDENA2F(N,A,B)
Input: N,A,B Output: A,B Sorts the real array A(N) into ascending numerical order. The additional array B(N) is simultaneously changed in parallel with the array A(N). Note that both input arrays are returned rearranged. We follow the Heapsort method, as described by D. Knuth, The Art of Computer Programming (pag.146, 5.2.3.). INTEGER N -> input number of data in A REAL A(N) -> data matrix to be sorted REAL B(N) -> data matrix to be sorted in parallel with matrix A
SUBROUTINE ORDENA2I(N,A,B)
Input: N,A,B Output: A,B Sorts the integer array A(N) into ascending numerical order. The additional array B(N) is simultaneously changed in parallel with the array A(N). Note that both input arrays are returned rearranged. We follow the Heapsort method, as described by D. Knuth, The Art of Computer Programming (pag.146, 5.2.3.). INTEGER N -> input number of data in A INTEGER A(N) -> data matrix to be sorted INTEGER B(N) -> data matrix to be sorted in parallel with matrix A
SUBROUTINE POLFIT(X,Y,SIGMAY,NPTS,NTERMS,MODE,A,CHISQR)
>>> This subroutine is based on the subroutine from Data Reduction and
Error Analysis for the Physical Sciences (Bevington, 1969) <<<
LEAST-SQUARES FIT TO A POLYNOMIAL
INPUT: X - ARRAY FOR INDEPENDENT VARIABLE
Y - ARRAY FOR DEPENDENT VARIABLE
SIGMAY - STANDARD DEVIATIONS FOR Y DATA POINTS
NPTS - NUMBER OF PAIRS OF DATA POINTS
NTERMS - NUMBER OF COEFFICIENTS (DEGREE + 1)
MODE - METHOD OF WEIGHTING (0 = NO WEIGHTING)
OUTPUT:A - ARRAY OF COEFFICIENTS
CHISQR - REDUCED CHI SQUARE FOR FIT
IT USES FUNCTION DETERM TO EVALUATE DETERMINANT OF MATRIX
SUPPORTS NTERM UP TO 20
FOR DETAILS SEE BEVINGTON(1969)
SUBROUTINE PSEUDOFIT(XF,YF,NF,NTERMS,YRMSTOL,WEIGHT,POWER,LUP,A)
Input: XF,YF,NF,NTERMS,YRMSTOL,WEIGHT,POWER,LUP
Output: A
Calculate the polynomial fit to the upper/lower side of a set of data
points.
REAL XF(NF),YF(NF) -> data points to be fitted
INTEGER NF -> number of data points
INTEGER NTERMS -> number of coeffcients
REAL YRMSTOL -> stopping criterion for DOWNHILL
REAL WEIGHT -> weighting factor to enhance one side of the fit
REAL POWER -> power to be used to compute distances
LOGICAL LUP -> .TRUE.: fit upper side
.FALSE.: fit lower side
REAL A(NTERMS) -> fitted coefficients
REAL FUNCTION RANRED(NSEED)
Input: NSEED
Output: RANRED (function)
Return a random number in the range [0,1) using the intrinsic fortran
function RAND(). If NSEED<0 a previous call to SRAND(TIME()) is also
performed.
INTEGER NSEED -> NSEED=0: RANRED returns the next random number in the
sequence.
NSEED<0: RANRED performs a previous call to the
intrinsic fortran function SRAND(TIME()), and
generates a random number in the new sequence.
In this case NSEED returns 0.
NSEED>0: RANRED performs a previous call to the
intrinsic fortran function SRAND(NSEED), and
generates a random number in the new sequence.
In this case NSEED returns 0.
SUBROUTINE REBINING(X,Y,N,XX,YY,M,XINI,XINC)
Input: X,Y,N,XINI,XINC Output: XX,YY,M Rebin a data table X(1:N),Y(1:N). If X(J)-X(J-1) .GT. XINC the routine performs a linear interpolation. REAL X(N) -> ordered input array (not necesarilly equally-spaced) REAL Y(N) -> input array INTEGER N -> no. of points in input array REAL XX(M) -> equally-spaced output array REAL YY(M) -> output array INTEGER M -> no. of points in output array REAL XINI -> equal to XX(1) REAL XINC -> equal to XX(J)-XX(J-1) for all J
REAL FUNCTION RINTERP(LAMBDA)
Input: LAMBDA Output: RINTERP (function) Calculate the extinction value A(lambda)/E(B-V) for a given wavelength LAMBDA. The tabulated data correspond to Savage & Mathis (1979, Ann. Rev. Astron. Astrophys., 17, 13). REAL LAMBDA -> input wavelength
SUBROUTINE RVREBIN(RADVEL,NCHAN,S,SS,STWV,DISP)
Input: RADVEL,NCHAN,S,STWV,DISP Output: SS This subroutine applies a radial velocity shift to a spectrum. REAL RADVEL -> radial velocity to be applied INTEGER NCHAN -> number of channels REAL S(NCHAN) -> initial spectrum REAL SS(NCHAN) -> shifted spectrum REAL STWV -> central wavelength of the first pixel REAL DISP -> dispersion (Angs./pixel) in the wavelength direction
SUBROUTINE CUBSPL__(X,Y,N,IMODE,S,A,B,C)
IMPLICIT NONE
INTEGER N
REAL X(N),Y(N)
INTEGER IMODE
REAL S(N)
REAL A(N),B(N),C(N)
local parameters
INTEGER NMAX
PARAMETER (NMAX=100)
local variables
INTEGER I,I1,I2
REAL DX1,DX2,DY1,DY2
REAL SS(0:NMAX,4)
REAL H
-----------------------------------------------------------------------------
verificamos que al menos tenemos dos puntos
IF(N.LT.2)THEN
WRITE(*,100)'FATAL ERROR in subroutine CUBSPL__: '
WRITE(*,100)' no. of data points too small: '
WRITE(*,*)N
STOP
END IF
verificamos que el numero de puntos no sea demasiado grande
IF(N.GT.NMAX)THEN !en realidad vale con N-1, pero lo dejamos asi
WRITE(*,100)'FATAL ERROR in subroutine CUBSPL__: '
WRITE(*,100)' no. of data points too large: '
WRITE(*,*)N
STOP
END IF
verificamos que IMODE toma un valor posible
IF((IMODE.LT.1).OR.(IMODE.GT.4))THEN
WRITE(*,100)'FATAL ERROR in subroutine CUBSPL__: '
WRITE(*,100)' invalid IMODE value:'
WRITE(*,*)IMODE
STOP
END IF
-----------------------------------------------------------------------------
calculamos la matriz reducida, que contiene N-2 filas x 4 columnas
DX1=X(2)-X(1)
DY1=6.*(Y(2)-Y(1))/DX1
DO I=1,N-2
DX2=X(I+2)-X(I+1)
DY2=6.*(Y(I+2)-Y(I+1))/DX2
SS(I,1)=DX1
SS(I,2)=2.*(DX1+DX2)
SS(I,3)=DX2
SS(I,4)=DY2-DY1
DX1=DX2
DY1=DY2
END DO
-----------------------------------------------------------------------------
modificamos la primera y ultima fila segun el valor de IMODE
IF(IMODE.EQ.1)THEN !no hay que modificar nada en este caso
.............................................................................
ELSEIF(IMODE.EQ.2)THEN
SS(1,2)=SS(1,2)+(X(2)-X(1))
SS(N-2,2)=SS(N-2,2)+(X(N)-X(N-1))
.............................................................................
ELSEIF(IMODE.EQ.3)THEN
DX1=X(2)-X(1)
DX2=X(3)-X(2)
SS(1,2)=(DX1+DX2)*(DX1+2.*DX2)/DX2
SS(1,3)=(DX2*DX2-DX1*DX1)/DX2
DX1=X(N-1)-X(N-2)
DX2=X(N)-X(N-1)
SS(N-2,1)=(DX1*DX1-DX2*DX2)/DX1
SS(N-2,2)=(DX2+DX1)*(DX2+2.*DX1)/DX1
.............................................................................
ELSEIF(IMODE.EQ.4)THEN
notar que en este caso tambien cambia el taman~o de la matriz
DX1=X(2)-X(1)
DY1=(Y(2)-Y(1))/DX1
SS(0,1)=1.
SS(0,2)=2.*DX1
SS(0,3)=DX1
SS(0,4)=6.*(DY1-S(1))
DX1=X(N)-X(N-1)
DY1=(Y(N)-Y(N-1))/DX1
SS(N-1,1)=DX1
SS(N-1,2)=2.*DX1
SS(N-1,3)=0.
SS(N-1,4)=6.*(S(N)-DY1)
.............................................................................
END IF
-----------------------------------------------------------------------------
dimensiones de la matriz a resolver
IF(IMODE.EQ.4)THEN
I1=1
I2=N-1
ELSE
I1=2
I2=N-2
END IF
resolvemos el sistema tridiagonal, eliminando en primer lugar los elementos
que se encuentran por debajo de la diagonal
DO I=I1,I2
SS(I,1)=SS(I,1)/SS(I-1,2)
SS(I,2)=SS(I,2)-SS(I,1)*SS(I-1,3)
SS(I,4)=SS(I,4)-SS(I,1)*SS(I-1,4)
END DO
y ahora hacemos la sustitucion desde atras
SS(I2,4)=SS(I2,4)/SS(I2,2)
DO I=I2-1,I1-1,-1
SS(I,4)=(SS(I,4)-SS(I,3)*SS(I+1,4))/SS(I,2)
END DO
-----------------------------------------------------------------------------
los valores de las derivadas segundas se almacenan en la matriz S
DO I=I1-1,I2
S(I+1)=SS(I,4)
END DO
IF(IMODE.EQ.1)THEN
S(1)=0.
S(N)=0.
ELSEIF(IMODE.EQ.2)THEN
S(1)=S(2)
S(N)=S(N-1)
ELSEIF(IMODE.EQ.3)THEN
DX1=X(2)-X(1)
DX2=X(3)-X(2)
S(1)=((DX1+DX2)*S(2)-DX1*S(3))/DX2
DX1=X(N-1)-X(N-2)
DX2=X(N)-X(N-1)
S(N)=((DX1+DX2)*S(N-1)-DX2*S(N-2))/DX1
ELSEIF(IMODE.EQ.4)THEN !no hay nada que cambiar
END IF
-----------------------------------------------------------------------------
finalmente calculamos los coeficientes
DO I=1,N-1
H=X(I+1)-X(I)
A(I)=(S(I+1)-S(I))/(6.*H)
B(I)=S(I)/2.
C(I)=(Y(I+1)-Y(I))/H-H*(2.*S(I)+S(I+1))/6.
END DO
-----------------------------------------------------------------------------
0 FORMAT(A,$)
END
*****************************************************************************
SUBROUTINE CUBSPLX__(X,Y,A,B,C,N,I0,X0,Y0)
Input: X,Y,A,B,C,N,I0,X0
Output: Y0
The subroutine returns the cubic spline evaluated at X0, using the
coefficients determined in a previous call to CUBSPL__. The spline defined in
the interval between X(I),Y(I) and X(I+1),Y(I+1) is given by:
Y = A(I)*(X-X(I))**3 + B(I)*(X-X(I))**2 + C(I)*(X-X(I)) + D(I)
If X0.LT.X(1), I=1 is employed (first computed spline)
If X0.GT.X(N), I=N-1 is employed (last computed spline)
REAL X(N) -> X-values fitted with CUBSPL__
REAL Y(N) -> Y-values fitted with CUBSPL__
REAL A(N) -> spline coefficients
REAL B(N) -> spline coefficients
REAL C(N) -> spline coefficients
INTEGER N -> number of data points
INTEGER I0 -> initial location to start the search of the place of X0 in
the X array
REAL X0 -> X-value where the spline function will be evaluated
REAL Y0 -> spline value at X0
SUBROUTINE SELBANDS(CBAND,NPBAND,WV,RES)
Input: CBAND Output: NPBAND,WV,RES Return the response curves of some common photometric bands. CHARACTER*1 CPBAND -> photometric band: U,B,V INTEGER NPBAND -> no. of points which define the output table REAL WV(NPBANDMAX) -> wavelengths REAL RES(NPBANDMAX) -> response curve
SUBROUTINE SELINDEX(NINDEX,WV,FWV,ITI,CLABEL)
Input: NINDEX
Output: WV,FWV,ITI,CLABEL
Return the bandpass limits of atomic, and molecular indices (and the D4000).
The subroutine looks first for a file called 'myindex.dat' in the current
directory. If this file does not exist, the program then looks for a file
called 'index.dat' (located in the subdirectory 'files' of the distribution
package). If this last file is also missing, the program stops.
INTEGER NINDEX -> index number. If NINDEX=0 the routine returns ITI
with the total number of defined indices.
REAL WV(NWVMAX) -> wavelength limits.
REAL FWV(NWVMAX/4) -> constant factors to be applied to the data in
the absorption bands.
INTEGER ITI -> index type:
ITI = -??: slope
ITI = 1 : molecular
= 2 : atomic
= 3 : D4000
= 4 : B4000
= 5 : color
= ????: generic with
ITI= C x 100 + L, where C=No. continuum regions
L=No. absorption regions
Cmin=1, Cmax=99, Lmin=1, Lmax=99
(ITImin=101, ITImax=9999)
CHARACTER*8 CLABEL -> character string with index identification
SUBROUTINE SELLINES(NTYPE,NLINES,WV,CLABEL)
Input: NTYPE
Output: NLINES,WV,CLABEL
Return the wavelength location of typical lines.
INTEGER NTYPE -> type of lines:
NTYPE=0: Balmer serie
NTYPE=1: typical emission lines
NTYPE=2: typical sky lines
NTYPE=3: typical absorption lines
INTEGER NLINES -> no. of returned lines
REAL WV(NLINMAX) -> wavelengths
CHARACTER*8 CLABEL(NLINMAX) -> character string with line identification
SUBROUTINE SHINDEX(LINDOK,MODE)
Input: LINDOK,MODE
Show a list of available indices (with an 80 character width format) in the
subroutine SELINDEX.
LOGICAL LINDOK(NINDMAX) -> if .TRUE. the index is shown in the list
INTEGER MODE -> if MODE=0 two extra options are displayed, namely
-1:EXIT and 0:ALL.
SUBROUTINE SPLFIT(N,X,Y,ND,XD,YRMSTOL,NOUT,XOUT,YOUT,XMIN,XMAX,SIGMA,LPLOTS)
Input: N,X,Y,ND,XD,YRMSTOL,NOUT,XMIN,XMAX,SIGMA,LPLOTS Output: XOUT,YOUT Least-squares fit to splines, using ND knots located at XD(). Input data are X(N), Y(N). XOUT(NOUT), YOUT(NOUT) are the output values which are computed in the range from XMIN to XMAX. The knot location determines the X(),Y() range employed in the fit (which is performed in the interval from XD(1) to XD(ND)) INTEGER N -> initial number of points in input data REAL X(N) -> sorted input data REAL Y(N) -> input data INTEGER ND -> number of knots REAL XD(ND) -> X location of each knot REAL YRMSTOL -> stopping criterion for DOWNHILL INTEGER NOUT -> number of points in output REAL XOUT(NOUT) -> output data REAL YOUT(NOUT) -> output data REAL XMIN -> = XOUT(1) REAL XMAX -> = XOUT(NOUT) REAL SIGMA -> sigma of the fit LOGICAL LPLOTS -> if .TRUE. some plots are performed
SUBROUTINE SUBPRECE(TII,RAI,DECI,TFF,RAF,DECF)
Input: TII,RAI,DECI,TFF Output: RAF,DECF Transformation of coordinates given for an equinox to another equinox (precession effect). DOUBLE PRECISION TII -> initial equinox (year) DOUBLE PRECISION RAI -> initial right ascension (hours) DOUBLE PRECISION DECI -> initial declination (degrees) DOUBLE PRECISION TIF -> final equinox (year) DOUBLE PRECISION RAF -> final right ascension (hours) DOUBLE PRECISION DECF -> final declination (degrees)
SUBROUTINE ULOGREB(COPC,S,N,CRVAL,CRPIX,CDELT,SS,M,STWV,DISP)
Input: COPC,S,N,CRVAL,CRPIX,CDELT,STWV,DISP
Output: SS,M
Transform a spectrum S(N) in logarithmic wavelength scale into a linear
wavelength scale.
CHARACTER*1 COPC -> type of wavelength calibration of input spectrum
COPC='1': CTYPE='WAVE' (linear)
COPC='2': CTYPE='WAVE-LOG' (log10)
COPC='3': CTYPE='WAVE-LOG' (ln)
COPC='4': CTYPE='wavenumber'
REAL S(N) -> input spectrum
INTEGER N -> no. of points in input spectrum
REAL CRVAL,CRPIX,CDELT -> wavelength calibration of input spectrum
if COPC='4', CRVAL and CRPIX are wavenumbers
REAL SS(M) -> output spectrum
INTEGER M -> no. of points in output spectrum
REAL STWV, DISP -> linear wavelength calibration for output spectrum
|
|