/* @(#)moscalib.c 17.1.1.1 (ESO-IPG) 01/25/02 17:54:47 */ /* @(#)moscalib.c 1.6 (ESO_IPG) 2/17/94 14:49:34 */ /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++ */ /* .COPYRIGHT (C) 1993 European Southern Observatory */ /* .IDENT lncalib.c */ /* .AUTHORS Pascal Ballester (ESO/Garching) */ /* Cristian Levin (ESO/La Silla) */ /* Sabine Moehler (LSW) */ /* .KEYWORDS Spectroscopy, Long-Slit, MOS */ /* .PURPOSE Wavelength Calibration for MOS Spectra */ /* .VERSION 1.0 Package Creation 17-MAR-1993 */ /* 1.1 modified for MOS 20-AUG-1993 */ /* */ /* .INPUT/OUTPUT: IN_A : line table (def. linpos.tbl) IN_B : line catalog MOS : mos table (slitlet positions) (def. mos) OUT_A : output coefficients table (def. coerbr.tbl) INPUTC : mode (Ident, Linear or Fors)(Constant/Variable fit) INPUTD/D/1/6 : start(1),start(2),step(1),step(2),startw,pixel INPUTR/R/1/3 : alpha,maxdev,tol INPUTI/I/1/5 : degx,degl,miniter,maxiter,ystart INPUTI(10) : Debug level ------------------------------------------------------- */ #include #include #include #include /*#include */ #include #include #include int inull; /* Representation of NULL values */ float rnull; double dnull; char mode[2]; int recall; #define DEGY 3 #define DEGXY 2 /*---------------------declaration of prototypes-----------------------------*/ #ifdef __STDC__ int fit_select(double [], double [], double [], int, double, int [], double [], double [], int, double [], int, int); double comp_dif (double [], double [], double [], int); double mode_init(char, double [], double [], double [], int, int); int auto_id( int, int, int, double [], int[], float [], int [], double [], int [], int [], double, int [], int [], int *, double [], int); int write_icol(int, int, int [], int, int []); int read_select( int, int, int []); int read_ident(double *, double *, int, double *, double *, int *); int fit_select_2D(double [], double [], double [], double [], int, double, int [], double [], double [], double [], int, double [], int); int auto_2d( int, int, int, double [], int[], float [], int [], double [], int [], int [], double, int [], int [], int *, double [], int); int read_ident_2D(double *, double *, double *, int, double *, double *, double *, int *); #else int fit_select(); double comp_dif (); double mode_init(); int auto_id(); int write_icol(); int read_select(); int read_ident(); int fit_select_2D(); int auto_2d(); int read_ident_2D(); #endif main() { char line[80], lincat[60], outtab[60]; /*mode[2]; */ char mod_save[1], mos[60], text[120]; /*Table names*/ int tidlin, tidcat, tidmos; /* Table id numbers */ int nbcol, nbrow, nsort, allcol, allrow; /* Descriptors of line */ int nbcolcat, nbrowcat, nsortcat, allcolcat, allrowcat; /* Descriptors of lincat */ int nbcolmos, nbrowmos, nsortmos, allcolmos, allrowmos; int catwav, linwav, linwavc, linslit, mosoff, mosslit; int linx, liny, lindif, numlin[8], numcat[3]; /* Column numbers */ int linrej, *reject, *select, *yrow, nyrow, nbsel; int ipar[10], degx, minter, maxter, ystart; int maxslit, prevslit=0, *slitrow, numslit, islit, islit1; int actvals, kunit, null; int i, j, k, slit_id, ystart_row, rowstart, rowend, direction; int slitsel, iref, itest, disp; int slit, y_index, row, shift_tol; int tmps, ix, rownum, rownum1, idrow, nid, loop=0; float rpar[10], alpha, maxdev, tol; /* Real parameters */ double *colx, *off; double mos_save[100], *coly, *yval, *slitval, *colslit; double delta, mindelta, spar[11], pixel; double xpos[50], ident[50], dpar[11], start[2], step[2]; double chi, *xid, *yid, *lid, xide[50], lide[50], Xid[50], Lid[50]; double tmpo, xtmp[100], diffx, ytmp; SCSPRO("moscalib"); /* Read name of table line in keyword IN_A, lincat in IN_B and mode in INPUTC(1:2). Mode can be Ident, Linear, or Fors and Constant or Variable linelist*/ SCKGETC("IN_A", 1, 60, &actvals, line); SCKGETC("IN_B", 1, 60, &actvals, lincat); SCKGETC("IN_C", 1, 60, &actvals, mos); SCKGETC("OUT_A", 1, 60, &actvals, outtab); SCKGETC("INPUTC",1, 2, &actvals, mode); SCKRDR("INPUTR", 1, 3, &actvals, rpar, &kunit, &null); SCKRDI("INPUTI", 1, 5, &actvals, ipar, &kunit, &null); SCKRDI("INPUTI", 10, 1, &actvals, &disp,&kunit, &null); /* Display level */ SCKRDD("INPUTD", 1, 7, &actvals, dpar, &kunit, &null); mod_save[0] = toupper(mode[0]); if (toupper(mode[0]) == 'I') { SCKRDD("XPOS", 1, 50, &actvals, xpos, &kunit, &null); SCKRDD("LID", 1, 50, &actvals, ident, &kunit, &null); } degx = ipar[0]; minter = ipar[2]; maxter = ipar[3]; ystart = ipar[4]; alpha = rpar[0]; maxdev = rpar[1]; tol = rpar[2]; start[0] = dpar[0]; start[1] = dpar[1]; step[0] = dpar[2]; step[1] = dpar[3]; dpar[10] = 0.; /* linear dispersion */ /* Open table mos and search for column :xoffset. All table are mapped to double precision arrays */ if (TCTOPN(mos, F_I_MODE, &tidmos)) SCTPUT("**** Error while opening table mos"); TCIGET(tidmos, &nbcolmos, &nbrowmos, &nsortmos, &allcolmos, &allrowmos); numslit = nbrowmos; /* Finds column :SLIT in table mos */ TCCSER(tidmos, ":SLIT", &mosslit); if (mosslit == (-1)) SCTPUT("**** Column :SLIT not found"); /* Finds column :XOFFSET in table mos */ TCCSER(tidmos, ":XOFFSET", &mosoff); if (mosoff == (-1)) SCTPUT("**** Column :XOFFSET not found"); /* Open table lincat and search for column :WAVE. All table are mapped to double precision arrays */ if (TCTOPN(lincat, F_I_MODE, &tidcat)) SCTPUT("**** Error while opening line catalog"); TCIGET(tidcat, &nbcolcat, &nbrowcat, &nsortcat, &allcolcat, &allrowcat); /* Finds column :WAVE in line catalog */ TCCSER(tidcat, ":WAVE", &catwav); if (catwav == (-1)) SCTPUT("**** Column :WAVE not found"); /* Storing adresses in array numcat[3] */ numcat[0] = tidcat; numcat[1] = catwav; numcat[2] = nbrowcat; /* Open table line and search for columns :X, :Y, :WAVE, :WAVEC, :RESIDUAL, :REJECT. All table columns are mapped to double precision arrays */ if (TCTOPN(line, F_IO_MODE, &tidlin)) SCTPUT("**** Error while opening table line.tbl"); TCIGET(tidlin, &nbcol, &nbrow, &nsort, &allcol, &allrow); TCCSER(tidlin, ":X", &linx); if (linx == (-1)) SCTPUT("**** Column :X not found"); TCCSER(tidlin, ":Y", &liny); if (liny == (-1)) SCTPUT("**** Column :Y not found"); TCCSER(tidlin, ":SLIT", &linslit); if (linslit == (-1)) SCTPUT("**** Column :SLIT not found"); TCCSER(tidlin, ":WAVE", &linwav); if (linwav == (-1)) TCCINI(tidlin, D_R8_FORMAT, 1, "F10.3", "Angstrom", "WAVE", &linwav); TCCSER(tidlin, ":WAVEC", &linwavc); if (linwavc == (-1)) TCCINI(tidlin, D_R8_FORMAT, 1, "F10.3", "Angstrom", "WAVEC", &linwavc); TCCSER(tidlin, ":RESIDUAL", &lindif); if (lindif == (-1)) TCCINI(tidlin, D_R8_FORMAT, 1, "F10.3", "Angstrom", "RESIDUAL", &lindif); TCCSER(tidlin, ":REJECT", &linrej); if (linrej == (-1)) TCCINI(tidlin, D_I4_FORMAT, 1, "I6", "Rejection Code", "REJECT", &linrej); /* Storing adresses in array numlin[8] */ numlin[0] = tidlin; numlin[1] = linwav; numlin[2] = linwavc; numlin[3] = lindif; numlin[4] = linrej; numlin[5] = linx; numlin[6] = liny; numlin[7] = nbrow; /* Read other parameters: INPUTI(1) Degree in X INPUTI(2) Degree in Y INPUTI(3) Minimal number of iterations INPUTI(4) Maximal number of iterations INPUTR(1) Parameter Alpha (tolerance) for matching (allowed residual for line identification = alpha*(minimum residual)(between the matched and the next catalog line OR between the searched and the next computed line)) INPUTR(2) Maximal standard deviation of residuals (pixels) INPUTR(3) Tolerance of individual lines for fitting (pixels if >0, wavel. units if <0 */ if (disp >= 50) { SCTPUT(text); sprintf(text,"Line table : %s ",line); SCTPUT(text); sprintf(text,"Line catalog : %s ",lincat); SCTPUT(text); sprintf(text,"Mode : %s ",mode); SCTPUT(text); sprintf(text,"Nb of iterations (min,max) : %d , %d",minter,maxter); SCTPUT(text); sprintf(text,"Degree : %d ",degx); SCTPUT(text); sprintf(text,"Tolerance for individual lines (pixels) : %f",tol); SCTPUT(text); sprintf(text,"Rejection parameter Alpha : %f",alpha); SCTPUT(text); sprintf(text,"Maximum mean deviation (pixels) : %f",maxdev); SCTPUT(text); } /* Initialisations */ /* Get representation of NULL values and allocate memory for auxiliary arrays */ TCMNUL(&inull, &rnull, &dnull); /* Load line table and line catalog in memory. Table columns are stored in double although they can be of any type (I, R*4 or R*8) */ select = (int *) osmmget((nbrow+1)*sizeof(int)); nbsel = read_select(tidlin, nbrow, select); sprintf(text," Number of lines (total, selected) : %d, %d", nbrow, nbsel); SCTPUT(text); colx = (double *) osmmget((nbrow+1)*sizeof(double)); read_col(tidlin, nbrow, linx, colx, dnull); coly = (double *) osmmget((nbrow+1)*sizeof(double)); read_col(tidlin, nbrow, liny, coly, dnull); colslit = (double *) osmmget((nbrow+1)*sizeof(double)); read_col(tidlin, nbrow, linslit, colslit, dnull); /* find out how many slitlets are selected */ slitsel = 0; iref = -1; for (i=1; i<=nbsel; i++) { itest = colslit[i]; if (itest != iref) { slitsel++; iref = colslit[i]; } } maxslit = slitsel; sprintf(text," number of slitlets : %i",maxslit); SCTPUT(text); off = (double *) osmmget((nbrowmos+1)*sizeof(double)); /* For table MOS the selection flag is ignored */ /* Get offsets as function of slit number */ for (i = 1; i <= nbrowmos; i++) off[i] = 0; for (i = 1; i <= nbrowmos; i++) { TCERDD(tidmos, i, mosoff, &tmpo, &null); TCERDI(tidmos, i, mosslit, &tmps, &null); off[tmps] = tmpo; } slitval = (double *) osmmget((nbrow+1)*sizeof(double)); yval = (double *) osmmget((nbrow+1)*sizeof(double)); reject = (int *) osmmget((nbrow+1)*sizeof(int)); slitrow = (int *) osmmget((nbrowmos+2)*sizeof(int)); yrow = (int *) osmmget((nbrow+1)*sizeof(int)); xid = (double *) osmmget((nbrow+1)*sizeof(double)); yid = (double *) osmmget((nbrow+1)*sizeof(double)); lid = (double *) osmmget((nbrow+1)*sizeof(double)); /* Go through the line table to identify the starting row number for the different y values */ y_index = 1, row=0; while (y_index <= nbsel) { /* Aendern: yval direkt aus tidlin einlesen und dann nur wenn yval = coly yrow beschreiben. so muesset Einhalten der Selektion moeglich sein */ yval[++row] = coly[y_index]; yrow[row] = y_index; while (y_index <= nbsel && coly[y_index] == yval[row]) ++y_index; } nyrow = row; yrow[nyrow+1] = nbsel+1; /* Go through the line table to identify the starting row number for the different slitlets */ row = 1; slit = 0; while (row <= nyrow) { y_index = yrow[row]; slitval[++slit] = colslit[y_index]; iref = colslit[y_index]; slitrow[iref] = row; while (y_index <= nbsel && colslit[y_index] == slitval[slit]) { ++row; y_index = yrow[row]; } } slitval[slitsel+1] = iref+1; /* y-position of the selected slitlets - - TS*/ slitrow[iref+1] = nyrow + 1; /* coerbr-table: row number in the table to write the coefs! - TS*/ /* Note: if the coefs are calculated for slitlet 1 in six rows, the data of the second slitlet starts in row slitrow[2] = 7 array counts from 1 to (number of slitlets + 1) - TS*/ /* For mode IDENT: Find the position ystart_row corresponding to the position of ystart For all other modes: Set ystart = first row of first selected slitlet */ if (toupper(mode[0]) == 'I') { /* ystart is in world coordiantes */ for (row=1; row<=nyrow; row++) { delta = yval[row] - ystart; if (delta < 0) delta = delta*(-1.); if (row == 1) mindelta = delta; if (delta <= mindelta) { mindelta = delta; ystart_row = row; y_index = yrow[row]; slit_id = colslit[y_index]; off[0] = off[slit_id]; } } } else { ystart = coly[1]; for (row=1; row<=nyrow; row++) { delta = yval[row] - ystart; if (delta < 0) { delta = delta*(-1.); } if (row == 1) mindelta = delta; if (delta <= mindelta) { mindelta = delta; ystart_row = row; y_index = yrow[row]; slit_id = colslit[y_index]; off[0] = off[slit_id]; dpar[6] = off[slit_id]; } } } /* Estimate a first dispersion relation in mode ident or read it in modes linear and fors */ recall = 0; if (toupper(mode[1]) == 'T') setrefdeg_2D(degx,DEGY,DEGXY); /*DEGY=4 DEGXY=1 from the header*/ pixel = mode_init(mode[0], xpos, ident, dpar, degx, nbrow); if (pixel == -1) goto ciao; if (toupper(mode[0]) == 'R') recall--; /* save display value */ ipar[9] = disp; /* identify the selected slitlet only */ for (islit = 1; islit <= maxslit; islit++) /*loop all slitlets to select*/ { ipar[5] = slitval[islit]; if (slitval[islit] == slit_id) /*select here*/ { slitsel--; for (direction=1; direction>=-1; direction-=2) /*two directions*/ { if (direction == 1) { setrefdeg(degx); mos_initdisp(outtab,"NEW",1); /* islit1: ID-number of next slitlet -TS0896*/ /* rowstart, rowend: 1st and last row in a slitlet */ rowstart = ystart_row; islit1 = slitval[islit+1]; rowend = slitrow[islit1]-1; if (toupper(mode[1]) == 'V') { auto_id(rowstart, rowend, direction, yval, yrow, rpar, ipar, dpar, numlin, numcat, pixel, select, reject, &loop, mos_save, nyrow); finishdisp(); } else { auto_2d(rowstart, rowend, direction, yval, yrow, rpar, ipar, dpar, numlin, numcat, pixel, select, reject, &loop, mos_save, nyrow); setrefdeg(degx); finishdisp(); } } if (direction == -1 && ystart_row > 1) { setrefdeg(degx); mos_initdisp(outtab,"OLD",1); rowstart = ystart_row; islit1 = slitval[islit]; rowend = slitrow[islit1]; if (toupper(mode[1]) == 'V') { auto_id(rowstart, rowend, direction, yval, yrow, rpar, ipar, dpar, numlin, numcat, pixel, select, reject, &loop, mos_save, nyrow); finishdisp(); } else { /*here I need the mode_init to start with 1D-Coefs*/ pixel = mode_init(mode[0],xpos,ident,dpar,degx,nbrow); auto_2d(rowstart, rowend, direction, yval, yrow, rpar, ipar,dpar, numlin, numcat, pixel, select, reject, &loop, mos_save, nyrow); setrefdeg(degx); finishdisp(); } } } /* End of loop on direction */ } /*end of (slitval[islit] == slit_id)*/ } /*end of for all slitlets*/ if (mos_save[1] > 0 && toupper(mode[0]) == 'R') { recall = 1; /* successful fit -- mode 'R' is active now */ /* mos_savedisp(mos_save);*//* save the first dispersion coeficients*/ } /* initialize spar[10] */ /* For mode LINEAR the initial parameters dpar[10] are kept */ for (i = 1; i <= 10; i++) spar[i] = dpar[i]; for (i = 1; i <= maxslit; i++) /* Loop for all selected slitlets */ { /* i <= maxslit */ islit = (int) slitval[i]; ipar[5] = islit; if ( (islit != slit_id && slitsel > 0) || toupper(mode[0]) == 'I' ) { /* islit != slit_id && slitsel > 0 -- slit_id has been done before*/ slitsel--; spar[6] = off[islit]; /* Estimate a first dispersion relation */ if (toupper(mode[0]) == 'L') { /*linear no recall - the coefs are reinitialized*/ setdisp(degx,mos_save); pixel = mode_init(mode[0], xpos, ident, spar, degx, nbrow); if (pixel == -1) goto ciao; } if (toupper(mode[0]) == 'I' || mod_save[0] == 'I') /*interactive*/ { /*interactive mode*/ nid = 0; /* read identifications and correct for offset */ for (j=0; j<50; j++) { if (ident[j] != 0) { Xid[nid] = xpos[j]-off[slit_id]+off[islit]; Lid[nid++] = ident[j]; } } /* read x-positions of lines in first row of resp. slitlet */ ix=0; idrow = slitrow[islit]; rownum = yrow[idrow]; rownum1 = yrow[idrow+1]; while (rownum < rownum1) xtmp[ix++] = colx[rownum++]; /* identify lines */ SCKRDI ("SHIFTTOL", 1, 1, &actvals, &shift_tol, &kunit, &null); for (j = 0; j < 50; j++) lide[j] = 0; for (j = 0; j < 50; j++) xide[j] = 0; j = 0; k = 0; while (j < nid) { for (k = 0; k < ix; k++) { diffx = xtmp[k]-Xid[j]; if (diffx < 0) diffx = -diffx; if (diffx <= shift_tol) { lide[k] = Lid[j]; xide[k] = xtmp[k]; } } j++; } pixel = mode_init(mode[0], xide, lide, spar, degx, nbrow); } /*end of interactive mode*/ direction = 1; setrefdeg(degx); mos_initdisp(outtab,"OLD",1); rowstart = slitrow[islit]; islit1 = slitval[i+1]; rowend = slitrow[islit1]-1; ipar[0] = degx; if (toupper(mode[1]) != 'V') { auto_2d(rowstart, rowend, direction, yval, yrow, rpar, ipar, spar, numlin, numcat, pixel, select, reject, &loop, mos_save, nyrow); setrefdeg(degx); finishdisp(); } else { auto_id(rowstart, rowend, direction, yval, yrow, rpar, ipar, spar, numlin, numcat, pixel, select, reject, &loop, mos_save, nyrow); finishdisp(); } prevslit = islit; } /*end of islit != slitid*/ else prevslit = slitval[i-1]; } /*end of i<=maxslit*/ ciao: TCSINI(tidcat); TCTCLO(tidcat); TCSINI(tidlin); TCTCLO(tidlin); osmmfree((char *) off); osmmfree((char *)slitval ); osmmfree((char *)yval ); osmmfree((char *)reject ); osmmfree((char *)slitrow ); osmmfree((char *)yrow ); osmmfree((char *)select ); osmmfree((char *)colx ); osmmfree((char *)coly ); osmmfree((char *)colslit ); osmmfree((char *) xid); osmmfree((char *) yid); osmmfree((char *) lid); SCSEPI(); } /*--------------------------------------------------------------------------*/ #ifdef __STDC__ double mode_init(char mod1, double x[], double id[], double linpar[], int degx, int nbrow) #else double mode_init( mod1, x, id, linpar, degx, nbrow ) char mod1; double x[],id[],linpar[]; int degx, nbrow; #endif /* Estimate a first dispersion relation in mode ident or read it in modes linear and fors */ { char text[120]; int i, nid; double pixel, *xid, *lid, chi; xid = (double *) osmmget((nbrow+1)*sizeof(double)); lid = (double *) osmmget((nbrow+1)*sizeof(double)); switch (toupper(mod1)) { case ('I') : { nid = 0; for (i=0; i<50; i++) { if (id[i] != 0) { xid[++nid] = x[i]; lid[nid] = id[i]; } } if (nid < 2) { sprintf (text,"Not enough identifications... Exiting.\n"); SCTPUT(text); goto ciao_m; } set_zero(degx); pixel = mos_fit_disp (&nid, °x, xid, lid, &chi); osmmfree((char *) xid); osmmfree((char *) lid); return (pixel); break; } case ('L') : { double coefs[2]; coefs[0] = linpar[4]-linpar[6]*linpar[5];/* x = lam0-offs*disp */ coefs[1] = linpar[5];/*disp*/ /*for (i=2; i<=6; i++) coefs[i] = 0;*//* these are degrees in Y*/ setdisp(1,coefs);/* degree=1, coef[i] = coefs[i-1]*/ pixel = linpar[5]; osmmfree((char *) xid); osmmfree((char *) lid); return (pixel); break; } case ('R') : { double coefs[2]; coefs[0] = linpar[4]-linpar[6]*linpar[5];/* x = lam0-offs*disp */ coefs[1] = linpar[5];/*disp*/ /*for (i=2; i<=6; i++) coefs[i] = 0;*/ setdisp(1,coefs); /* degree=1, coef[i] = coefs[i-1]*/ pixel = linpar[5]; osmmfree((char *) xid); osmmfree((char *) lid); return (pixel); break; } default: { osmmfree((char *) xid); osmmfree((char *) lid); (void) sprintf(text, "Error in moscalib.c: Unknown calibration method %c\n",mod1); SCETER(9,text); /* already leaves this function... */ break; } } ciao_m: pixel = -1.; osmmfree((char *) xid); osmmfree((char *) lid); return (pixel); } /*---------------------------------------------------------------------------*/ #ifdef __STDC__ auto_id( int rstart, int rend, int direc, double val_y[], int row_y[], float rlist[], int ilist[], double dlist[], int nlin[], int ncat[], double pixel, int sel[], int rej[], int *loop, double mos_save[], int nyrow) #else auto_id( rstart, rend, direc, val_y, row_y, rlist, ilist, dlist, nlin, ncat, pixel, sel, rej, loop, mos_save, nyrow) int rstart,rend,direc,row_y[],ilist[],nlin[],ncat[],sel[],rej[],*loop,nyrow; float rlist[]; double val_y[],pixel,mos_save[],dlist[]; #endif { char text[120]; int tidlin, linwavc, lindif, linrej, linwav, linx, liny, nbrow; int *select, *sline, itmax, itmin, tidcat, nbrowcat, catwav, dgx; int nblines, nbsel, y_index, iter, prevmatch, end_of_iter, slit; int ungraceful_exit, linb, nmatch, verif=0, cal=0; int lindeg = 1, i, nid, row, kunit, disp; double *colcat, *coldif, *colx, *coly, *colwav, *colwavc; double *xline, *yline, *xid, *lid, *xtmp, start[2], step[2]; double slit_save[100], chi, current_y, stdres, stdpix, tolwav; float alpha, maxdev, tol; #ifdef __STDC__ int match(int, double [], double [], double [], double [], int, double [], int, double, double *, double, int []); #else int match(); #endif /* Read parameters */ dgx = ilist[0]; itmin = ilist[1]; itmax = ilist[2]; disp = ilist[9]; start[1] = dlist[1]; step[1] = dlist[3]; alpha = rlist[0]; maxdev = rlist[1]; tol = rlist[2]; /* Get representation of NULL values and allocate memory for auxiliary arrays */ TCMNUL(&inull, &rnull, &dnull); /* read stored table adresses: nlin = numlin -- ncat = numcat*/ tidlin = nlin[0]; linwav = nlin[1]; linwavc = nlin[2]; lindif = nlin[3]; linrej = nlin[4]; linx = nlin[5]; liny = nlin[6]; nbrow = nlin[7]; tidcat = ncat[0]; catwav = ncat[1]; nbrowcat = ncat[2]; /* Load line table and line catalog in memory. Table columns are stored in double although they can be of any type (I, R*4 or R*8) */ select = (int *) osmmget((nbrow+1)*sizeof(int)); nbsel = read_select(tidlin, nbrow, select); colcat = (double *) osmmget((nbrowcat+1)*sizeof(double)); read_col(tidcat, nbrowcat, catwav, colcat, dnull); colx = (double *) osmmget((nbsel+1)*sizeof(double)); read_col(tidlin, nbrow, linx, colx, dnull); coly = (double *) osmmget((nbsel+1)*sizeof(double)); read_col(tidlin, nbrow, liny, coly, dnull); coldif = (double *) osmmget((nbsel+1)*sizeof(double)); read_col(tidlin, nbrow, lindif, coldif, dnull); xline = (double *) osmmget((nbsel+1)*sizeof(double)); yline = (double *) osmmget((nbsel+1)*sizeof(double)); sline = (int *) osmmget((nbsel+1)*sizeof(int)); colwav = (double *) osmmget((nbsel+1)*sizeof(double)); colwavc = (double *) osmmget((nbsel+1)*sizeof(double)); xid = (double *) osmmget((nbsel+1)*sizeof(double)); lid = (double *) osmmget((nbsel+1)*sizeof(double)); rej = (int *) osmmget((nbsel+1)*sizeof(int)); xtmp = (double *) osmmget((nbsel+1)*sizeof(double)); /*xtmp=xline-off*/ /* Write null values in all output columns */ if (loop == 0) { for (row=0; row<=nbrow; row++) { sline[row] = row; colwav[row] = dnull; colwavc[row] = dnull; coldif[row] = dnull; rej[row] = inull; } write_dcol(tidlin, nbrow, sline, linwav, colwav); write_dcol(tidlin, nbrow, sline, linwavc,colwavc); write_dcol(tidlin, nbrow, sline, lindif, coldif); write_icol(tidlin, nbrow, sline, linrej, rej); *loop=1; } /* Read a set of values at first row */ current_y = val_y[rstart]; /* Read slitlet number */ slit = ilist[5]; nblines = 0; for (y_index=row_y[rstart]; y_index= 50) { sprintf(text,"Variable dispersion for slit nr. %d ystart = %7.1f", ilist[5], current_y); SCTPUT(text); } while (!end_of_iter && !ungraceful_exit) { /* iterative match loop */ /* Compute wavelength for all lines: */ if (recall != 0 && iter == 0) /*recall disp... if (mode[0] = 'R')*/ { if (recall == -1) /* first linear fit for first slitlet*/ { mos_eval_disp(xline, colwavc, nblines); } else /* recall disp... for other slitlets*/ { setdisp(dgx,mos_save); mos_eval_disp(xtmp, colwavc, nblines); } } else mos_eval_disp(xline, colwavc, nblines); /* match (nblines) lines against the catalog: the matched lines of the first row are recalled in the loop*/ nmatch = match(verif, colwav, colwavc, yline, coldif, nblines, colcat, nbrowcat, alpha, &stdres, dnull, rej); stdpix = stdres/pixel; if (disp >= 100) { sprintf(text," row Y = %4d: matching %2d lines out of %2d", (int)current_y, nmatch, nblines); SCTPUT(text); } iter++; /* Test end of iteration: */ end_of_iter = ((iter >= itmax) || (iter > itmin && nmatch == prevmatch)); prevmatch = nmatch; if (!(ungraceful_exit = (stdpix > maxdev))) /*not crash condition*/ { /* Read the matched lines -- if !ungraceful_exit */ read_ident(xline, colwav, nblines, xid, lid, &nid); pixel = mos_fit_disp (&nid, &dgx, xid, lid, &chi); if (pixel < 0.) ungraceful_exit = TRUE; } } /* end_of_iter */ /* Fit the single rows... of LINPOS-table*/ if (ungraceful_exit) { /*ungraceful_exit*/ sprintf(text,"Sorry, wrong identifications...\n"); SCTPUT(text); set_zero(dgx); stdres = -1.; for (row=rstart; row <= rend+direc ; row+=direc) { linb = (int) ((current_y - start[1])/step[1] + 1.5); mos_writedisp (row, -1, linb, current_y, nyrow, stdres); cal = -1; SCKWRI("CAL", &cal, slit, 1, &kunit); } if (recall != 0) { /*reset coefs to the start values*/ setdisp(1,mos_save); setrefdeg(dgx); } } else { /*not ungraceful_exit*/ mos_savedisp(slit_save); /* ... slitsave[i] = coef[i+1]*/ cal = 1; SCKWRI("CAL", &cal, slit, 1, &kunit); /* Fit dispersion coefficients for all rows: */ /* changed from != to <= for the case rstart = rend -- TS*/ for (row=rstart; row != rend+direc ; row+=direc) { /*loop on the rows*/ /* Read a set of values at constant y */ current_y = val_y[row]; nblines = 0; for (y_index=row_y[row]; y_index0.) ? tol*pixel : -tol; nid = fit_select (xline, colwav, coldif, nblines, tolwav, rej, xid, lid, nid, colwavc, dgx, (int) current_y); if (nid >= 2) pixel = mos_fit_disp (&nid, &dgx, xid, lid, &chi); if (pixel > 0 && nid >= 2) { mos_eval_disp(xline, colwavc, nblines); stdres = comp_dif(colwav, colwavc, coldif, nblines); if (disp >= 100) { sprintf(text," Final selection for Y = %4d: %2d lines out of %2d", (int)current_y, nid, nblines); SCTPUT(text); } if (disp >= 50) { sprintf(text," RMS = %6.2f - Tolerance = %6.2f (wav. units)", stdres, tolwav); SCTPUT(text); if (disp < 100) disp = 40; } linb = (int) ((current_y - start[1])/step[1] + 1.5); mos_writedisp (row, slit, linb, current_y, nyrow, stdres); } if (pixel <= 0 || nid <= 1) { stdres = -1.; set_zero(dgx); mos_writedisp (row, -1, linb, current_y, nyrow, stdres); } write_dcol(tidlin, nblines, sline, linwav, colwav); write_dcol(tidlin, nblines, sline, linwavc,colwavc); write_dcol(tidlin, nblines, sline, lindif, coldif); write_icol(tidlin, nblines, sline, linrej, rej); if (pixel > 0 && nid > dgx && recall == -1) { mos_savedisp(slit_save); mos_savedisp(mos_save); recall = recall + 2; if (disp >= 50) { printf(" save mos_disp : "); printdisp(); } } } /* End of loop on row */ /* Perform linear fit to get a mean linear dispersion for the transformation*/ if (pixel > 0 && nid >= dgx && recall == 0) { mos_fit_disp (&nid, &lindeg, xid, lid, &chi); mos_savedisp(mos_save); /* ...mos_save[i] = coef[i+1]*/ } setrefdeg(dgx); } /* end of (un)graceful_exit*/ /* mos_savedisp(mos_save); *//* ...mos_save[i] = coef[i+1]*/ osmmfree((char *) select); osmmfree((char *) colcat); osmmfree((char *) colx); osmmfree((char *) coly); osmmfree((char *) coldif); osmmfree((char *) xline); osmmfree((char *) yline ); osmmfree((char *) sline ); osmmfree((char *) colwav ); osmmfree((char *) colwavc ); osmmfree((char *) xid ); osmmfree((char *) lid ); osmmfree((char *) rej); osmmfree((char *) xtmp); }/* end of AUTO_ID */ /*---------------------------------------------------------------------------*/ #ifdef __STDC__ auto_2d( int rstart, int rend, int direc, double val_y[], int row_y[], float rlist[], int ilist[], double dlist[], int nlin[], int ncat[], double pixel, int sel[], int rej[], int *loop, double mos_save[], int nyrow) #else auto_2d( rstart, rend, direc, val_y, row_y, rlist, ilist, dlist, nlin, ncat, pixel, sel, rej, loop, mos_save, nyrow) int rstart,rend,direc,row_y[],ilist[],nlin[],ncat[],sel[],rej[],*loop,nyrow; float rlist[]; double val_y[],pixel,mos_save[],dlist[]; #endif { char text[120]; int tidlin, linwavc, lindif, linrej, linwav, linx, liny, nbrow; int *select, *sline, itmax, itmin, tidcat, nbrowcat, catwav, dgx; int nblines, nbsel, y_index, iter, prevmatch, end_of_iter, slit; int ungraceful_exit, linb, nmatch, verif=0, cal=0; int lindeg = 1, i, nid, row, kunit, disp,ypixel[9999],numypix; double *colcat, *coldif, *colx, *coly, *colwav, *colwavc; double *xline, *yline, *xid, *yid, *lid, *xtmp, *ytmp, start[2], step[2]; double slit_save[100], chi, current_y, stdres, stdpix, tolwav; float alpha, maxdev, tol, offset; #ifdef __STDC__ int match(int, double [], double [], double [], double [], int, double [], int, double, double *, double, int []); #else int match(); #endif /* Read parameters */ dgx = ilist[0]; itmin = ilist[1]; itmax = ilist[2]; disp = ilist[9]; start[1] = dlist[1]; step[1] = dlist[3]; offset = dlist[6]; alpha = rlist[0]; maxdev = rlist[1]; tol = rlist[2]; /* Get representation of NULL values and allocate memory for auxiliary arrays */ TCMNUL(&inull, &rnull, &dnull); /* read stored table adresses : nlin=numlin ncat=numcat -- from main{}*/ tidlin = nlin[0]; linwav = nlin[1]; linwavc = nlin[2]; lindif = nlin[3]; linrej = nlin[4]; linx = nlin[5]; liny = nlin[6]; nbrow = nlin[7]; tidcat = ncat[0]; catwav = ncat[1]; nbrowcat = ncat[2]; /* Load line table and line catalog in memory. Table columns are stored in double although they can be of any type (I, R*4 or R*8) */ select = (int *) osmmget((nbrow+1)*sizeof(int)); nbsel = read_select(tidlin, nbrow, select); colcat = (double *) osmmget((nbrowcat+1)*sizeof(double)); read_col(tidcat, nbrowcat, catwav, colcat, dnull); colx = (double *) osmmget((nbsel+1)*sizeof(double)); read_col(tidlin, nbrow, linx, colx, dnull); coly = (double *) osmmget((nbsel+1)*sizeof(double)); read_col(tidlin, nbrow, liny, coly, dnull); coldif = (double *) osmmget((nbsel+1)*sizeof(double)); read_col(tidlin, nbrow, lindif, coldif, dnull); xline = (double *) osmmget((nbsel+1)*sizeof(double)); yline = (double *) osmmget((nbsel+1)*sizeof(double)); sline = (int *) osmmget((nbsel+1)*sizeof(int)); colwav = (double *) osmmget((nbsel+1)*sizeof(double)); colwavc = (double *) osmmget((nbsel+1)*sizeof(double)); xid = (double *) osmmget((nbsel+1)*sizeof(double)); yid = (double *) osmmget((nbsel+1)*sizeof(double)); lid = (double *) osmmget((nbsel+1)*sizeof(double)); rej = (int *) osmmget((nbsel+1)*sizeof(int)); xtmp = (double *) osmmget((nbsel+1)*sizeof(double)); ytmp = (double *) osmmget((nbsel+1)*sizeof(double)); /* Write null values in all output columns */ if (loop == 0) { for (row=0; row<=nbrow; row++) { sline[row] = row; colwav[row] = dnull; colwavc[row] = dnull; coldif[row] = dnull; rej[row] = inull; } write_dcol(tidlin, nbrow, sline, linwav, colwav); write_dcol(tidlin, nbrow, sline, linwavc,colwavc); write_dcol(tidlin, nbrow, sline, lindif, coldif); write_icol(tidlin, nbrow, sline, linrej, rej); *loop=1; } current_y = val_y[rstart]; /* Read a set of values at first row */ slit = ilist[5]; /* Read slitlet number */ nblines = 0; for (y_index=row_y[rstart]; y_index 5) { /* the real 2D-fit: */ if (disp >= 50 && iter == 0) { sprintf(text,"2D-Dispersion for slit nr. %d ystart = %7.1f numypix = %6d", ilist[5], current_y, numypix); SCTPUT(text); } if (toupper(mode[0]) == 'I' && iter == 0) { mos_savedisp(mos_save); /*read 1-D coefs*/ for (i=dgx+1;i<=dgx+6;i++) mos_save[i] = 0.; setdisp_2D(dgx,mos_save); /*write 2-D coefs*/ } /* Compute wavelength for all lines: */ if (recall != 0 && iter == 0) /*recall disp... if (mode[0] = 'F')*/ { if (recall == -1) /* first linear fit for first slitlet*/ { mos_eval_disp_2D(xtmp, yline, colwavc, nblines); } else /* recall disp... for other slitlets*/ { setdisp_2D(dgx,mos_save); mos_eval_disp_2D(xtmp, ytmp, colwavc, nblines); } } else mos_eval_disp_2D(xline, yline, colwavc, nblines); /* match all lines against the catalog */ nmatch = match(verif, colwav, colwavc, yline, coldif, nblines, colcat, nbrowcat, alpha, &stdres, dnull, rej); /* for (i=1;i<=nblines;i++) printf("xline %6.1f yline[i] %7.1f colwavc %8.3f colwav %8.3f \n",xline[i],yline[i],colwavc[i],colwav[i]);*/ stdpix = stdres/pixel; if (disp >= 100) { sprintf(text," row Y = %4d: matching %4d lines out of %4d", (int)current_y, nmatch, nblines); SCTPUT(text); } /* Test end of iteration: */ iter++; end_of_iter = ((iter >= itmax) || (iter > itmin && nmatch == prevmatch)); prevmatch = nmatch; if (!(ungraceful_exit = (stdpix > maxdev))) { /* Estimate the new dispersion relation on the matched lines */ read_ident_2D(xline, yline, colwav, nblines, xid, yid, lid,&nid); pixel = mos_fit_disp_2D (&nid, &dgx, xid, yid, lid, &chi); if (pixel < 0.) ungraceful_exit = TRUE; } } else { /* constant dispersion over the slitlet */ if (disp >= 50 && iter == 0) { sprintf(text,"constant-dispersion for slit nr. %d ystart = %7.1f numypix = %6d", ilist[5], current_y, numypix); SCTPUT(text); } /* Compute wavelength for all lines: */ if (recall != 0 && iter == 0) /*recall disp... if (mode[0] = 'F')*/ { if (recall == -1) /* first linear fit for first slitlet*/ { mos_eval_disp(xtmp, colwavc, nblines); } else /* recall disp... for other slitlets*/ { if (toupper(mode[1]) == 'T') { /* the coefs are still in 2D - but not if iter>0 */ setdisp_2D(dgx, mos_save); mos_eval_disp_2D(xtmp, ytmp, colwavc, nblines); } else { setdisp(dgx,mos_save); mos_eval_disp(xtmp, colwavc, nblines); } } } else mos_eval_disp(xline, colwavc, nblines); /* match all lines against the catalog: */ nmatch = match(verif, colwav, colwavc, yline, coldif, nblines, colcat, nbrowcat, alpha, &stdres, dnull, rej); stdpix = stdres/pixel; if (disp >= 100) { sprintf(text," row Y = %4d: matching %4d lines out of %4d", (int)current_y, nmatch, nblines); SCTPUT(text); } /* Test end of iteration: */ iter++; end_of_iter = ((iter>=itmax) || (iter>itmin && nmatch==prevmatch)); prevmatch = nmatch; if (!(ungraceful_exit = (stdpix > maxdev))) { /* Estimate the new dispersion relation on the matched lines */ read_ident(xline, colwav, nblines, xid, lid,&nid); pixel = mos_fit_disp(&nid, &dgx, xid, lid, &chi); if (pixel < 0.) ungraceful_exit = TRUE; } } } /* end_of_iter */ /* Fit the data in two dimensions... data from LINPOS-table*/ /* reject the lines with residuals greater then the limit*/ if (ungraceful_exit) { /*ungraceful_exit*/ sprintf(text,"\nSlitlet %3d: Sorry, wrong identifications...\n",ilist[5]); SCTPUT(text); set_zero(dgx); stdres = -1.; for (row=rstart; row != rend+direc ; row+=direc) { linb = (int) ((current_y - start[1])/step[1] + 1.5); mos_writedisp (row, -1, linb, val_y[row], nyrow, stdres); cal = -1; SCKWRI("CAL", &cal, slit, 1, &kunit); } if (recall == -1) { /*reset coefs to the start values*/ setdisp(1,mos_save); setrefdeg(dgx); } } else { /*(not un)graceful_exit*/ /* Fit dispersion coefficients for all rows: */ mos_savedisp(slit_save); /* ...slit_save[i] = coef[i+1]*/ cal = 1; SCKWRI("CAL", &cal, slit, 1, &kunit); current_y = val_y[row]; nblines = 0; /* loop on the rows -- Read values at all y for 2D-fit and constant: */ for (y_index=row_y[rstart]; y_index 5) { /* the real 2D-fit */ setdisp_2D(dgx,slit_save); /* coef[i] = slit_save[i-1]*/ mos_eval_disp_2D(xline, yline, colwavc, nblines); comp_dif(colwav, colwavc, coldif, nblines); tolwav = (tol>0.) ? tol*pixel : -tol; nid = fit_select_2D (xline, yline, colwav, coldif, nblines, tolwav, rej, xid, yid, lid, nid, colwavc, dgx); if (nid >= 2*numypix) pixel = mos_fit_disp_2D(&nid, &dgx, xid, yid, lid, &chi); if (pixel > 0 && nid >= 2*numypix) { /* here we have achieved convergence*/ mos_eval_disp_2D(xline, yline, colwavc, nblines); stdres = comp_dif(colwav, colwavc, coldif, nblines); if (disp >= 100) { sprintf(text," Final selection for Y = %4d: %2d lines out of %2d", (int)current_y, nid, nblines); SCTPUT(text); } if (disp >= 50) { sprintf(text," Slit = %3d RMS = %6.2f - Tolerance = %6.2f (wav. units)", ilist[5], stdres, tolwav); SCTPUT(text); if (disp < 100) disp = 40; } for (row = rstart; row != rend+direc ; row += direc) { /* save the coeficients*/ ypixel[row] = (int) ((val_y[row]-start[1])/step[1] + 1.5); mos_writedisp_2D(row,slit, ypixel[row],val_y[row],nyrow,stdres); } } if (pixel <= 0 || nid <= numypix) { /* here the real desaster */ stdres = -1.; set_zero(dgx); for (row = rstart; row != rend+direc ; row += direc) { /* save crash coeficients*/ ypixel[row] = (int) ((val_y[row]-start[1])/step[1] + 1.5); mos_writedisp_2D (row, -1, ypixel[row], val_y[row], nyrow, stdres); } setdisp(1,mos_save); setrefdeg(dgx); } write_dcol(tidlin, nblines, sline, linwav, colwav); write_dcol(tidlin, nblines, sline, linwavc,colwavc); write_dcol(tidlin, nblines, sline, lindif, coldif); write_icol(tidlin, nblines, sline, linrej, rej); /* Perform linear fit to get a mean linear dispersion for the transformation or set the initial parameters of the recall: */ if (pixel > 0 && nid > dgx*numypix && toupper(mode[0]) == 'L') { /* only for linear */ mos_fit_disp(&nid, &lindeg, xid, lid, &chi); mos_savedisp(mos_save); } /* method 2D-fit - set initial disp*/ if (pixel > 0 && nid >= dgx*numypix && recall == -1) { mos_savedisp(mos_save); recall = recall + 2; } setrefdeg_2D(dgx,DEGY,DEGXY); } /* end of the real 2D-fit*/ else { /* constant dispersion over the slitlet: */ setdisp(dgx,slit_save);/* coef[i] = slit_save[i-1]*/ mos_eval_disp(xline, colwavc, nblines); comp_dif(colwav, colwavc, coldif, nblines); tolwav = (tol>0.) ? tol*pixel : -tol; nid = fit_select(xline, colwav, coldif, nblines, tolwav, rej, xid, lid, nid, colwavc, dgx, (int) current_y); if (nid >= 2*numypix) pixel = mos_fit_disp(&nid, &dgx, xid, lid, &chi); if (pixel > 0 && nid >= 2*numypix) { mos_eval_disp(xline, colwavc, nblines); stdres = comp_dif(colwav, colwavc, coldif, nblines); if (disp >= 100 && iter == 0) { sprintf(text," Final selection for Y = %4d: %2d lines out of %2d", (int)current_y, nid, nblines); SCTPUT(text); } if (disp >= 50) { sprintf(text," Slit = %3d RMS = %6.2f - Tolerance = %6.2f (wav. units)", ilist[5], stdres, tolwav); SCTPUT(text); if (disp < 100) disp = 40; } for (row = rstart; row != rend+direc ; row += direc) { ypixel[row] = (int) ((val_y[row]-start[1])/step[1] + 1.5); mos_writedisp(row,slit,ypixel[row],val_y[row],nyrow,stdres); } } if (pixel <= 0 || nid <= numypix) { stdres = -1.; set_zero(dgx); for (row = rstart; row != rend+direc ; row += direc) { ypixel[row] = (int) ((val_y[row]-start[1])/step[1] + 1.5); mos_writedisp(row, -1, ypixel[row],val_y[row],nyrow,stdres); } } write_dcol(tidlin, nblines, sline, linwav, colwav); write_dcol(tidlin, nblines, sline, linwavc,colwavc); write_dcol(tidlin, nblines, sline, lindif, coldif); write_icol(tidlin, nblines, sline, linrej, rej); if (pixel > 0 && nid >= dgx*numypix && toupper(mode[0]) == 'L') { /* method linear */ mos_fit_disp (&nid, &lindeg, xid, lid, &chi); mos_savedisp(mos_save); } if (pixel > 0 && nid >= dgx*numypix && recall == -1) { /* method constant - set initial coefs*/ mos_savedisp(mos_save); recall = recall + 2; } setrefdeg(dgx); }/* end of method constant*/ } /* end of (un)graceful_exit*/ osmmfree((char *) select); osmmfree((char *) colcat); osmmfree((char *) colx); osmmfree((char *) coly); osmmfree((char *) coldif); osmmfree((char *) xline); osmmfree((char *) xtmp); osmmfree((char *) ytmp); osmmfree((char *) yline ); osmmfree((char *) sline ); osmmfree((char *) colwav ); osmmfree((char *) colwavc ); osmmfree((char *) xid ); osmmfree((char *) yid ); osmmfree((char *) lid ); osmmfree((char *) rej); }/* end of AUTO_2D */ /*---------------------------------------------------------------------------*/ #ifdef __STDC__ write_icol(int tid, int nb, int select[], int colnb, int col[]) #else write_icol( tid, nb, select, colnb, col ) int tid,nb,select[],colnb,col[]; #endif { int i; for (i=1; i<=nb; i++) TCEWRI(tid, select[i], colnb, &col[i]); } /*---------------------------------------------------------------------------*/ #ifdef __STDC__ double comp_dif ( double colwav[], double colwavc[], double coldif[], int nb) #else double comp_dif ( colwav, colwavc, coldif, nb) double colwav[],colwavc[],coldif[]; int nb; #endif { int i, number; double stdres; stdres = 0.; number = 0; for (i=1; i<=nb; i++) { if (colwav[i] != dnull) { coldif[i] = colwavc[i]-colwav[i]; number++; stdres += coldif[i]*coldif[i]; } } stdres = sqrt(stdres/(double) number); return stdres; } /*---------------------------------------------------------------------------*/ #ifdef __STDC__ int fit_select (double xline[], double colwav[], double coldif[], int nb, double tolwav, int reject[], double xid[], double lid[], int nid, double colwavc[], int dgx, int rowpos) #else int fit_select ( xline, colwav, coldif, nb, tolwav, reject, xid, lid, nid, colwavc, dgx, rowpos) double xline[],colwav[],coldif[],tolwav,xid[],lid[],colwavc[]; int nb,reject[],nid,dgx,rowpos; #endif { char text[120]; int maxpos, i=0; double px, maxres, resabs, chi, *tmpcolwav; /* double stdres; */ /* tmpcolwav construction to reject lines temporary row by row only note that colwav is used for any row of a slitlet and the match would be lost with any line rejected by fit_select -- TS0896*/ tmpcolwav = (double *) osmmget ((nb+1)*sizeof(double)); for (i=1; i<= nb; i++) tmpcolwav[i] = colwav[i]; maxres = tolwav; while (maxres >= tolwav) { maxres = 0.; nid = 0; for (i=1; i<=nb; i++) { if (reject[i] != RESIDUAL_GT_TOL && tmpcolwav[i] > 0) { resabs = (coldif[i] < 0) ? -coldif[i] : coldif[i]; if (resabs > maxres) maxpos = i, maxres = resabs; } } if (maxres > tolwav) { if ( tmpcolwav[maxpos] > 0.) { sprintf(text," bad line at %10.3f - residual: %8.3f (wav. units)", tmpcolwav[maxpos], maxres); SCTPUT(text); } tmpcolwav[maxpos] = dnull; reject[maxpos] = RESIDUAL_GT_TOL; read_ident(xline, tmpcolwav, nb, xid, lid, &nid); px = mos_fit_disp(&nid, &dgx, xid, lid, &chi); if (px > 0) { mos_eval_disp(xline, colwavc, nb); comp_dif(tmpcolwav, colwavc, coldif, nb); maxres=tolwav; } } else { for (i=1; i<=nb; i++) { if (reject[i] != RESIDUAL_GT_TOL && tmpcolwav[i] != dnull && xline[i]) { xid[++(nid)] = xline[i]; lid[nid] = tmpcolwav[i]; } } } } /* end of while */ osmmfree((char *) tmpcolwav); return nid; } /*---------------------------------------------------------------------------*/ #ifdef __STDC__ int fit_select_2D (double xline[],double yline[], double colwav[], double coldif[], int nb, double tolwav, int reject[],double xid[], double yid[], double lid[], int nid, double colwavc[], int dgx) #else int fit_select_2D ( xline, yline, colwav, coldif, nb, tolwav, reject, xid, yid, lid, nid, colwavc, dgx) double xline[],yline[],colwav[],coldif[],tolwav,xid[],yid[],lid[],colwavc[]; int nb,reject[],nid,dgx; #endif { char text[120]; int maxpos, i=0; double px, maxres, resabs, chi, lidrej; /* double stdres; */ maxres = tolwav; while (maxres >= tolwav) { maxres = 0.; nid = 0; lidrej = 0.; for (i=1; i<=nb; i++) { if (reject[i] != RESIDUAL_GT_TOL && colwav[i] != dnull) { resabs = (coldif[i] < 0) ? -coldif[i] : coldif[i]; if (resabs > maxres) maxpos = i, maxres = resabs; } } if (maxres > tolwav) { /*warning for identified lines only: */ sprintf(text," bad line at %10.3f - residual: %8.3f (wav. units)", colwav[maxpos], maxres); SCTPUT(text); colwav[maxpos] = dnull; reject[maxpos] = RESIDUAL_GT_TOL; read_ident_2D(xline, yline, colwav, nb, xid, yid, lid, &nid); px = mos_fit_disp_2D(&nid, &dgx, xid ,yid, lid, &chi); if (px > 0) { mos_eval_disp_2D(xline, yline, colwavc, nb); comp_dif(colwav, colwavc, coldif, nb); maxres = tolwav; } } else { for (i=1; i<=nb; i++) { if (reject[i] != RESIDUAL_GT_TOL && colwav[i] != dnull) { xid[++(nid)] = xline[i]; yid[nid] = yline[i]; lid[nid] = colwav[i]; } } } /* end of if */ } /* end of while */ return nid; } /* end of fit_select_2D */ /*---------------------------------------------------------------------------*/ #ifdef __STDC__ read_select( int tid, int nb, int col[]) #else read_select( tid, nb, col ) int tid,nb,col[]; #endif { int i, k=0, sel; for (i=1; i<=nb; i++) { TCSGET(tid, i, &sel); if (sel) col[++k] = i; } return(k); } /*---------------------------------------------------------------------------*/ #ifdef __STDC__ read_ident(double *colx, double *colid, int nbsel, double *xid, double *lid, int *nid) #else read_ident( colx, colid, nbsel, xid, lid, nid) double *colx,*colid,*xid,*lid; int nbsel,*nid; #endif { int i; *nid = 0; for (i=1; i<=nbsel; i++) { if (colid[i] != dnull && colx[i] !=dnull) { xid[++(*nid)] = colx[i]; lid[*nid] = colid[i]; } } }/*---------------------------------------------------------------------------*/ #ifdef __STDC__ read_ident_2D(double *colx, double *coly, double *colid, int nbsel, double *xid, double *yid, double *lid, int *nid) #else read_ident_2D( colx, coly, colid, nbsel, xid, yid, lid, nid) double *colx,*coly,*colid,*xid,*yid,*lid; int nbsel,*nid; #endif { int i; *nid = 0; for (i=1; i<=nbsel; i++) { if (colid[i] != dnull && colx[i] !=dnull) { xid[++(*nid)] = colx[i]; yid[*nid] = coly[i]; lid[*nid] = colid[i]; } } }