/* @(#)necoffset.c 17.1.1.1 (ES0-DMD) 01/25/02 17:51:30 */ /*=========================================================================== Copyright (C) 1995 European Southern Observatory (ESO) This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 675 Massachusetss Ave, Cambridge, MA 02139, USA. Corresponding concerning ESO-MIDAS should be addressed as follows: Internet e-mail: midas@eso.org Postal address: European Southern Observatory Data Management Division Karl-Schwarzschild-Strasse 2 D 85748 Garching bei Muenchen GERMANY ===========================================================================*/ /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++ */ /* .COPYRIGHT (C) 1995 European Southern Observatory */ /* .IDENT necoffset.c */ /* .AUTHOR Pascal Ballester (ESO/Garching) */ /* .KEYWORDS Spectroscopy, Echelle, Offset */ /* .PURPOSE Determines the slit offset of an exposure */ /* .VERSION 1.0 Creation 16-MAR-1995 */ /* */ /* .INPUT/OUTPUT: IN_A : Image name IN_B : Order table name INPUTC(11:10): Coefficients descriptor generic name INPUTI/I/1/4 : range, nb_positions, order nb min and max INPUTI/I/9/1 : Seed for random number generator INPUTI(10) : Debug level (0,1) ---------------------------------------------------------- */ #include #include #include #include #define ipos(col,row,siz) row * siz + col /* Pointer value */ #define nint(f) (int) (f + 0.5) /* Nearest integer */ #ifndef RAND_MAX #define RAND_MAX 32767 #endif double position(); main() { char frame[61], ordtab[61], coeff[11], descr[12]; int range, nb_positions, order_min, order_max, parami[4]; int dbx, nd, kmax, lmax; double coeffd[100], x, y; int coeffi[10], m; unsigned int seed; float *pntra, p_mean; int imnoa, naxis, npix[2], knul, stat, npar, kun; int iav, actvals, tid; int ncol, nrow, nsort, allcol, allrow, row; int ordcol, xcol, ycol, fitcol, rancol, xcent; double start[2], step[2]; char ident[TEXT_LEN+1], cunit[16*3+1]; int np, npts, sample, pos[2], x_pix, y_pix, p_max; int offs, pixel, *list; double incr, value; SCSPRO("offset"); SCKGETC ("IN_A", 1, 60, &actvals, frame); SCKGETC ("IN_B", 1, 60, &actvals, ordtab); SCKGETC ("INPUTC",11, 10,&actvals, coeff); SCKRDI ("INPUTI", 1, 4, &actvals, parami, &kun, &knul); range=parami[0], nb_positions=parami[1]; order_min=parami[2], order_max=parami[3]; SCKRDI ("INPUTI", 9, 2, &actvals, parami, &kun, &knul); seed = 2*parami[0]+1, dbx = parami[1]; srand(seed); /* Open Input image */ strcpy(ident, " "); strcpy(cunit, " "); SCIGET (frame, D_R4_FORMAT, F_I_MODE, F_IMA_TYPE, 2, &naxis, npix, start, step, ident, cunit, (char **)&pntra, &imnoa); nrow = npix[1]; ncol = npix[0]; /* Read Coefficients from order table */ TCTOPN (ordtab, F_IO_MODE, &tid); TCIGET (tid, &ncol, &nrow, &nsort, &allcol, &allrow); strcpy(descr,coeff); strcat(descr,"I"); SCDRDI(tid, descr, 4, 4, &iav, coeffi, &kun, &knul); kmax = coeffi[2], lmax = coeffi[3]; nd = (kmax+1)*(lmax+1); strcpy(descr,coeff); strcat(descr,"D"); SCDRDD(tid, descr, 1, nd, &iav, coeffd, &kun, &knul); np = (order_max - order_min +1)*npix[0]-1; npts = nb_positions; incr = (double)np/ (double)npts; list = (int *) osmmget((npts+1)*sizeof(int)); p_mean = 0.; for (sample=1; sample<=npts; sample++) { /* Pick random column and order number */ pick: pos[0] = rand_nb(order_min, order_max); pos[1] = rand_nb(1,npix[0]); x = start[0] + (pos[1]-1)*step[0]; y = position(x,pos[0],kmax,lmax,coeffd); y_pix = (y-start[1])/step[1]; x_pix = pos[1] - 1; if ((y_pix-range) < 0) goto pick; if ((y_pix+range) > (npix[1]-1)) goto pick; /* Determine offset of maximum signal */ for (offs=(-range); offs<=range; offs++) { pixel = ipos(x_pix, (y_pix+offs), npix[0]); if (offs == (-range)) {value = pntra[pixel]; p_max=offs;} else if (pntra[pixel] > value) {value = pntra[pixel]; p_max=offs;} } list[sample] = p_max; } sort(npts,list); pos[0] = nint(npts*0.375); pos[1] = nint(npts*0.625); p_mean = 0.; for (sample=pos[0]; sample<=pos[1]; sample++) p_mean += list[sample]; p_mean /= (double) (pos[1]-pos[0]+1.); /* printf("Found offset=%f\n",p_mean); */ SCKWRR("OUTPUTR", &p_mean, 1, 1, &kun); TCTCLO(tid); SCSEPI(); } int rand_nb(min,max) int min,max; { double frac; int result; frac = (double)rand()/(double)(RAND_MAX); if (frac > 1.) frac -= (int) frac; frac *= (max-min); result = min + nint(frac); return(result); } double position(x,m, kmax, lmax, dpar) double x; int m, kmax, lmax; double dpar[]; { int l, ip = -1, k; double y=0., dc=1., ww[100], yval=0.; for (l=0; l<=lmax; l++) { ip += 1, ww[ip] = dc, yval += ww[ip]*dpar[ip]; for (k=1; k<=kmax; k++) { ip += 1, ww[ip] = ww[ip-1]*x; yval += ww[ip]*dpar[ip]; } dc *= m; } return(yval); } /* Sorting algorithm: Heapsort method */ /* From: Numerical Recipes, Cambridge University Press, p. 226 */ sort(n,ra) int n; int ra[]; { int l,j,ir,i; float rra; l=(n >> 1)+1; ir=n; for (;;) { if (l > 1) rra=ra[--l]; else { rra=ra[ir]; ra[ir]=ra[1]; if (--ir == 1) { ra[1]=rra; return; } } i=l; j=l << 1; while (j <= ir) { if (j < ir && ra[j] < ra[j+1]) ++j; if (rra < ra[j]) { ra[i]=ra[j]; j += (i=j); } else j=ir+1; } ra[i]=rra; } }