! @(#)echregr.prg 3.3 (ESO-IPG) 4/8/92 18:10:06
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++
!
!.COPYRIGHT (C) 1991 European Southern Observatory
!.IDENT echregr.prg
!.AUTHOR Pascal Ballester, ESO - Garching
!.KEYWORDS Spectroscopy, Echelle,
!.PURPOSE Fit mesured position of orders by a bivariate polynomial
!.VERSION 1.0 Creation 13-OCT-1983
!
!-------------------------------------------------------
!
IF P1 .NE. "?" WRITE/KEYWORD DEFPOL/I/1/2 {P1}
DEFINE/PARAM P2 3 N "Number of outliers rejection:"
DEFINE/PARAM P3 2.0 N "Absolute tolerance (pixels):"
DEFINE/PARAM P4 4.5 N "Kappa-sigma factor :"
!
DEFINE/LOCAL CNT/I/1/1 0
DEFINE/LOCAL THRES/R/1/1 0.
DEFINE/LOCAL ORDER/I/1/1 0
DEFINE/LOCAL RMSMAX/R/1/1 0.
DEFINE/LOCAL RMS/R/1/1 0.
DEFINE/LOCAL REJECT/I/1/300 0 ALL
DEFINE/LOCAL NBREJ/I/1/1 0
!
WRITE/OUT "Warning: Present display shows the positions"
WRITE/OUT " used to fit the orders"
LOAD/ECHELLE ECHELLE RAW
!
SELECT/TABLE {ORDTAB} (:X.GE.1.AND.:X.LE.{IMSIZE(1)})
SELECT/TABLE {ORDTAB} SELECT.AND.(:Y.GE.{SCAN(1)}.AND.:Y.LE.{SCAN(2)})
WRITE/KEYWORD IN_A {ORDTAB}
WRITE/KEYWORD OUT_A middummw.tbl
INPUTI(1) = ECHORD(1)
RUN echregr
COPY/TABLE &w &s
SORT/TABLE &s :RMS
ORDER = ECHORD(1)/2
RMSMAX = 3.5*{middumms.tbl,:RMS,@{ORDER}}
WRITE/OUT "*** Order by order linear rms ***"
READ/TABLE &w
WRITE/OUT "Maximum admissible rms : {RMSMAX}"
DO ORDER = 1 {ECHORD(1)}
RMS = {middummw.tbl,:RMS,@{ORDER}}
IF RMS .GT. RMSMAX THEN
WRITE/OUT "Warning : Rejected order {ORDER}. Rms = {RMS}
NBREJ = NBREJ + 1
REJECT({NBREJ}) = ORDER
ENDIF
ENDDO
DO ORDER = 1 {NBREJ}
SELECT/TAB {ORDTAB} SELECT.AND.:ORDER.NE.{REJECT({ORDER})}
ENDDO
INPUTI(1) = DEFPOL(2) + 1 ! Degree Y + 1
INPUTI(2) = ECHORD(1) - NBREJ ! Number of orders minus rejected orders
IF INPUTI(1) .GT. INPUTI(2) THEN
WRITE/OUT "*****************************************************"
WRITE/OUT "**** Warning : Number of selected orders {INPUTI(2)}
WRITE/OUT "**** is too small for the current value of echelle
WRITE/OUT "**** parameter DEFPOL(2)={DEFPOL(2)}
WRITE/OUT "*****************************************************"\\
HELP/ECHELLE DEFPOL
ENDIF
REGRESSION/POLY {ORDTAB} :Y :X,:ORDER {DEFPOL(1)},{DEFPOL(2)} KEYLONG
SAVE/REGR {ORDTAB} COEFF KEYLONG
COMPUTE/REGR {ORDTAB} :YFIT = COEFF
COMPUTE/TABLE {ORDTAB} :RESIDUAL = :Y - :YFIT
THRES = {P4} * KEYLONGR(5)
IF THRES(1) .GT. {P3} THRES = {P3}
SELECT/TABLE {ORDTAB} ABS(:RESIDUAL).LT.{THRES}
IF {P2} .LT. 1 GOTO END
DO CNT = 1 {P2}
REGRESSION/POLY {ORDTAB} :Y :X,:ORDER {DEFPOL(1)},{DEFPOL(2)} KEYLONG
SAVE/REGR {ORDTAB} COEFF KEYLONG
COMPUTE/REGR {ORDTAB} :YFIT = COEFF
COMPUTE/TABLE {ORDTAB} :RESIDUAL = :Y - :YFIT
THRES = {P4} * KEYLONGR(5)
IF THRES(1) .GT. {P3} THRES = {P3}
SELECT/TABLE {ORDTAB} ABS(:RESIDUAL).LT.{THRES}
ENDDO
!
END:
-PURGE {ORDTAB}
To run the application execute the procedure via:
Midas 001> @@ echregr defpol ninter absres kappa
/* @(#)echregr.c 3.2 (ESO-IPG) 4/2/92 10:49:52 */
/* Program : echregr.c */
/* Author : P. Ballester - ESO Garching */
/* Date : 17.03.92 */
/* */
/* Purpose : */
/* Performs a linear fit on independent orders */
/* and writes downs the rms in a table */
/* */
/* Input : */
/* IN _A Name of order table */
/* INPUTI(1) Number of orders */
/* OUT_A Name of output table */
/* */
#include <math.h>
#include <tbldef.h>
#include <midas_def.h>
main()
{
int kunit;
char outtab[TEXT_LEN], ordtab[TEXT_LEN], text[TEXT_LEN];
int order, nb_order, present_order, row;
int rmscol, ordcol;
int order_col, xcol, ycol;
int ncol, nrow, nsort, allcol, allrow;
int select, null, status;
int iav, actvals, tid, tidord;
double x, y, cnt, sx, sy, sx2, sxy, sy2;
double det, a, b, rms;
SCSPRO("echregr");
SCKGETC ("IN_A", 1, 60, &actvals, ordtab);
SCKGETC ("OUT_A", 1, 60, &actvals, outtab);
SCKRDI ("INPUTI", 1, 1, &actvals, &nb_order, &kunit, &null);
TCTOPN (ordtab, F_I_MODE, &tidord);
TCIGET (tidord, &ncol, &nrow, &nsort, &allcol, &allrow);
TCCSER (tidord, "ORDER", &order_col);
TCCSER (tidord, "X", &xcol);
TCCSER (tidord, "Y", &ycol);
TCTINI (outtab, F_TRANS, F_IO_MODE, 3, 100, &tid);
TCCINI (tid, D_R4_FORMAT, 1, "I6", " ", "ORDER", &ordcol);
TCCINI (tid, D_R4_FORMAT, 1, "F12.4", " ", "RMS", &rmscol);
row = 1;
for (order=1; order<=nb_order; order++) {
cnt=0., sx=0., sy=0., sx2=0., sxy=0., sy2 = 0.;
present_order = order; /* Just to start the while () */
while (present_order == order) {
TCSGET(tidord, row, &select);
if (select) {
TCERDD(tidord, row, xcol, &x, &null);
TCERDD(tidord, row, ycol, &y, &null);
cnt += 1., sx += x, sy += y, sx2 += x*x, sy2 += y*y, sxy += x*y;
}
if (row >= nrow) break;
TCERDI(tidord, ++row, order_col, &present_order, &null);
}
if (cnt >= 3) {
det = cnt*sx2 - sx*sx;
a = (sy*sx2 - sx*sxy)/det;
b = (cnt*sxy - sx*sy)/det;
rms = sqrt((sy2 - a*a*cnt - 2.*b*a*sx - b*b*sx2)/cnt);
}
else rms = 99999.;
TCEWRI(tid, order, ordcol, &order);
TCEWRD(tid, order, rmscol, &rms);
}
/* Close everything and good bye. */
TCTCLO(tidord), TCTCLO(tid), SCSEPI();
}