C @(#)center.for 17.1.1.1 (ESO-DMD) 01/25/02 17:39:58 C=========================================================================== C Copyright (C) 1995 European Southern Observatory (ESO) C C This program is free software; you can redistribute it and/or C modify it under the terms of the GNU General Public License as C published by the Free Software Foundation; either version 2 of C the License, or (at your option) any later version. C C This program is distributed in the hope that it will be useful, C but WITHOUT ANY WARRANTY; without even the implied warranty of C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the C GNU General Public License for more details. C C You should have received a copy of the GNU General Public C License along with this program; if not, write to the Free C Software Foundation, Inc., 675 Massachusetts Ave, Cambridge, C MA 02139, USA. C C Correspondence concerning ESO-MIDAS should be addressed as follows: C Internet e-mail: midas@eso.org C Postal address: European Southern Observatory C Data Management Division C Karl-Schwarzschild-Strasse 2 C D 85748 Garching bei Muenchen C GERMANY C=========================================================================== C PROGRAM CENTER C C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C C.IDENTIFICATION: C program CENTER version 2.70 850322 C K. Banse ESO - Garching C 3.00 860528 3.10 860603 3.20 871124 3.30 880112 C 3.80 900207 3.90 900727 3.95 910319 4.00 920228 C C.KEYWORDS: C Gaussian profiles, intensity moments C C.PURPOSE: C Compute the central position of an object via Gaussian fit or C first moments of intensity C C.ALGORITHM: C See P.B. Stetson , 1979 Astron. J., 84, 1149 C C.INPUT/OUTPUT: C The following keywords are used: C IN_A/C/1/60 in_frame or in_frame,in_table C = CURSOR for cursor input C = + for no cursor C ACTION/C/1/2 (1) centering method, C G: Gaussian fitting C M: first moments of intensity C I: 2-dim Gaussian fitting + get angle of major axis C (2) invert_flag C Y: yes, we invert data (Gaussian valleys) C N: no, we have normal Gaussian peaks C INPUTC/C/1/30 specs for data storage C = :label for table or table,:label (a) C = descr/DESCR for descriptor (b) C = ? for just display of data C INPUTC/C/31/2 (a) identifier flag for table, C = I, append an identifier, else no. C (b) append flag for descriptor/table, C = A, append data to descriptor C = ?, start again at position 1 in descr C C OUTPUTR/R/1/11 last displayed results C OUTPUTI/I/1/1 no. of loops in program C C.VERSIONS C 3.00 clean up + use IDI interfaces C 3.30 added option for qualifier UGAUSS (Mauro Ghigo) C 3.50 START + STEP are now DOUBLE PRECISION C add zoom window + select flag for table stuff C 4.10 take care of 1-dim frames + tables (but there is still more C work to do with the display of values ...) C 4.20 add `invert_flag' for Gaussian valleys (versus peaks as usual) C C 010105 last modif C C-------------------------------------------------------------------------- C IMPLICIT NONE C INTEGER FELEM,IAV,IDFLAG,INFLAG,OUTFL,IDYES INTEGER KSTAT !retstat from centering: <0, 0,1,2,3, >3 INTEGER N,NN,NMAL,NLIM,NCOOS,INVFLG INTEGER NAXIS,TRUNAX,NROW,OUTCOL,NOCURS,FIRTIM INTEGER IMNO,STAT,STAT1,STAT2,LLABL INTEGER NPIX(3),NPIXT(3),DRAWY(3) INTEGER INBUF(5),NCOU,ALOW(3),BLOW(3),HIGH(3) INTEGER SPIX(4),XYALP(3) INTEGER ICUR1(2),ICUR2(2),NW(2),XU,XL,YU,YL INTEGER INCOLN(17),OUCOLN(17),LRVALS INTEGER TABNUL(12),TIDA,TIDB INTEGER COOS(4),OVCON,INPUT(3),CPIX(4),BOX(4),OLDBOX(4) INTEGER NAPP,ONEDIM,COOFF,FRPIX(2) INTEGER UNI(1),NULO,MAPROW,TIMNO INTEGER MADRID(1) INTEGER*8 PNTRA,PNTRC C CHARACTER FRAME*60 CHARACTER DESCR*15,APPFLG*2,CACT*8 CHARACTER CUNIT*48,IDENT*72,CBUF*82,ACTION*4,METHOD*4 CHARACTER IDF*8,OLDIDF*8,LABL*8 CHARACTER CSTAT*32,INFO*50 CHARACTER INTAB*60,OUTTAB*60,TABUW*16,TABUNI(17)*16 CHARACTER LABEL(17)*16,LABEL1(17)*16,LABEL2(17)*16,GLABEL*16 C DOUBLE PRECISION STEP(3),START(3) C REAL RBUF(17),PCUR1(6),PCUR2(6) REAL ICENT REAL RVALS(17),RINF(8),RV REAL XOUT,YOUT,XERR,YERR,XSIG,YSIG,XFWHM,YFWHM REAL XOFF,YOFF,XCUTS,ANGLE,ANGSIG,AXMAJ,AXMIN C DOUBLE PRECISION DIN(4),DOUT(4) !MaxDim = 4 in wrldco.c DOUBLE PRECISION DVALS(10) C LOGICAL SELFLG C INCLUDE 'MID_INCLUDE:ST_DEF.INC' INCLUDE 'MID_INCLUDE:IDIDEV.INC' INCLUDE 'MID_INCLUDE:IDIMEM.INC' C COMMON /VMR/ MADRID C DATA INFO /'switch cursor(s) on - next time we exit... '/ C DATA LABEL1 /'XSTART ','YSTART ','XEND ','YEND ', + 'XCEN ','YCEN ','XERR ','YERR ', + 'XSIG ','YSIG ', + 'XFWHM','YFWHM','ICENT', + 'STATUS ',' ',' ','IDENT '/ DATA LABEL2 /'XSTART ','YSTART ','XEND ','YEND ', + 'XCEN ','YCEN ','XERR ','YERR ', + 'XSIG ','YSIG ', + 'AX_MAJ ','AX_MIN ','ICENT', + 'STATUS ','ANGLE ','ANGLE_SIG ','IDENT '/ DATA GLABEL /'X_POSITION '/ DATA INCOLN /17*0/, OUCOLN /17*0/ DATA TABUW /'World Coords '/ DATA TABUNI /17*' '/ C DATA IDF /'ID '/, OLDIDF /' '/ DATA SPIX /4*1/, DESCR /' '/ DATA IDENT /' '/, CUNIT /' '/ DATA XYALP /3*0/, NPIX /3*1/, NPIXT /3*1/ DATA ALOW /3*1/, BLOW /3*1/, HIGH /3*1/ DATA CACT /'NNYY?C0 '/ DATA COOS /-1,-1,-1,-1/ DATA NCOOS /0/, IDFLAG /0/, ONEDIM /0/ DATA LRVALS /11/ DATA INTAB /' '/, OUTTAB /' '/ DATA DVALS /10*0.0/, DIN /4*1.0/ DATA START /3*0.0/, STEP /3*1.0/ C INCLUDE 'MID_INCLUDE:ST_DAT.INC' C C set up MIDAS environment + enable automatic error abort CALL STSPRO('CENTER') NCOU = 0 DRAWY(1) = -1 !set draw flag to start DRAWY(2) = 0 DRAWY(3) = 0 FIRTIM = 0 TIDB = -1 !default to NO output table C C get method CALL STKRDC('ACTION',1,1,2,IAV,ACTION,UNI,NULO,STAT) IF (ACTION(1:1).EQ.'M') THEN METHOD(1:) = 'MOM ' ELSE IF (ACTION(1:1).EQ.'I') THEN METHOD(1:) = 'IQE ' ELSE METHOD(1:) = 'GAU ' ENDIF IF ((ACTION(2:2).NE.'Y') .AND. (ACTION(2:2).NE.'y')) THEN INVFLG = 0 !usual Gaussian with max at center ELSE INVFLG = 1 ENDIF C IF (METHOD(1:1).EQ.'I') THEN DO 21, IAV=1,17 LABEL(IAV)(1:) = LABEL2(IAV)(1:) 21 CONTINUE ELSE DO 25, IAV=1,17 LABEL(IAV)(1:) = LABEL1(IAV)(1:) 25 CONTINUE ENDIF C C get input options CALL STKRDC('P1',1,1,80,IAV,CBUF,UNI,NULO,STAT) NN = INDEX(CBUF,',') C IF (NN.GT.1) THEN !image,Option FRAME(1:) = CBUF(1:NN-1)//' ' INTAB(1:) = CBUF(NN+1:)//' ' NN = INDEX(INTAB,',') IF (NN.GT.1) THEN !Option = xfpix,yfpix CALL GENCNV(INTAB,1,2,FRPIX,RBUF,STEP,IAV) IF (IAV .LT. 2) THEN CALL STETER(3,'Invalid frame pixels given...') ENDIF INFLAG = 4 !INFLAG = 4 for frame pixels ELSE !Option = table CALL CLNTAB(INTAB,INTAB,0) CALL TBTOPN(INTAB,F_IO_MODE,TIDA,STAT) INFLAG = 2 !INFLAG = 2 for table ENDIF ELSE IF ((CBUF(1:6).EQ.'CURSOR') .OR. + (CBUF(1:6).EQ.'cursor')) THEN CALL STKRDC('IN_A',1,1,60,IAV,FRAME,UNI,NULO,STAT) INFLAG = 1 !INFLAG = 1 for cursor COOFF = 0 !cursor_off flag ELSE FRAME(1:) = CBUF(1:60) INFLAG = 3 !INFLAG = 3 for complete image CALL STTPUT + ('We work on crude center of input frame. ',STAT) ENDIF ENDIF C C open input frame + map it CALL CLNFRA(FRAME,FRAME,0) CALL STFOPN(FRAME,D_R4_FORMAT,0,F_IMA_TYPE,IMNO,STAT) CALL STDRDI(IMNO,'NAXIS',1,1,IAV,NAXIS,UNI,NULO,STAT) CALL STDRDI(IMNO,'NPIX',1,3,IAV,NPIX,UNI,NULO,STAT) TRUNAX = NAXIS !get true NAXIS IF (NAXIS.GT.2) THEN IF (NPIX(3).EQ.1) THEN IF (NPIX(2).EQ.1) THEN TRUNAX = 1 ELSE TRUNAX = 2 ENDIF ELSE CALL STTPUT('only first plane of input frame used...',STAT) TRUNAX = 2 ENDIF ELSE IF (NAXIS.EQ.2) THEN IF (NPIX(2).EQ.1) TRUNAX = 1 ENDIF IF (TRUNAX.EQ.1) ONEDIM = 1 C CALL STDRDD(IMNO,'START',1,2,IAV,START,UNI,NULO,STAT) CALL STDRDD(IMNO,'STEP',1,2,IAV,STEP,UNI,NULO,STAT) IF (INVFLG.EQ.1) + CALL STDRDR(IMNO,'LHCUTS',4,1,IAV,XCUTS,UNI,NULO,STAT) C C init frame pixls <-> world coords conversion CALL FPXWCO(0,IMNO,DIN,DOUT,STAT) IF (STAT .GT. 0) + CALL STETER(11,'problems with WCS of input frame...') C C map temporary work frames, start with max. 220 rows MAPROW = 220 IF (ONEDIM.EQ.1) THEN N = NPIX(1) ELSE N = NPIX(1)*MAPROW ENDIF CALL STFYMP(N,D_R4_FORMAT,TIMNO,PNTRA,STAT) N = 200 * 200 CALL STFXMP(N,D_R4_FORMAT,PNTRC,STAT) C C get out_specs + output_options CALL STKRDC('INPUTC',1,1,32,IAV,CBUF,UNI,NULO,STAT) APPFLG(1:2) = CBUF(31:32) CBUF(31:32) = ' ' CALL UPCAS(APPFLG,APPFLG) C IF (CBUF(1:1).EQ.'+') THEN OUTFL = 0 GOTO 500 !neither table nor descriptor ENDIF C NN = INDEX(CBUF,',D') !test for ,D IF (NN.LE.1) NN = INDEX(CBUF,',d') C C handle descriptor storage IF (NN.GT.1) THEN !it's a descriptor DESCR = CBUF(1:NN-1) IF (APPFLG(1:1).EQ.'A') THEN !see, if data is to be appended FELEM = -1 ELSE FELEM = 1 ENDIF OUTFL = 1 !OUTFL = 1 for descriptor GOTO 500 !skip following table stuff ENDIF C C ++++ C handle table storage C ++++ C IF ( ((APPFLG(1:1).EQ.'I').OR.(APPFLG(2:2).EQ.'I')) + .AND. (INFLAG.EQ.1) ) IDFLAG = 1 CALL CLNTAB(CBUF(1:30),OUTTAB,0) DO 90, IAV=1,12 TABUNI(IAV)(1:) = TABUW(1:) 90 CONTINUE TABUNI(13)(1:) = 'Intensity ' TABUNI(14)(1:) = 'Center_Status ' IF (ACTION(1:1).EQ.'I') THEN OUTCOL = 16 TABUNI(15)(1:) = TABUW(1:) TABUNI(16)(1:) = TABUW(1:) ELSE OUTCOL = 14 ENDIF C C either same table involved C CALL GENEQF(OUTTAB,INTAB,IAV) IF (IAV.EQ.1) THEN !results go to same table TIDB = TIDA OUTFL = 2 !OUTFL = 2 for same table C DO 100 IAV=5,6 !see, if :XCENT,:YCENT already there CALL TBLSER(TIDA,LABEL(IAV),OUCOLN(IAV),STAT) IF (OUCOLN(IAV).LE.0) !No. We have to define it... + CALL TBCINI(TIDA,D_R8_FORMAT,1,'E24.15',TABUNI(IAV), + LABEL(IAV),OUCOLN(IAV),STAT) 100 CONTINUE C DO 120 IAV=7,OUTCOL !test, if columns already there CALL TBLSER(TIDA,LABEL(IAV),OUCOLN(IAV),STAT) IF (OUCOLN(IAV).LE.0) !No. We have to define it... + CALL TBCINI(TIDA,D_R4_FORMAT,1,'E12.5',TABUNI(IAV), + LABEL(IAV),OUCOLN(IAV),STAT) 120 CONTINUE C CALL TBLSER(TIDA,LABEL(17),OUCOLN(17),STAT) IF (OUCOLN(17).LE.0) !No. We have to define it... + CALL TBCINI(TIDA,D_C_FORMAT,8,'A8',TABUNI(17),LABEL(17), + OUCOLN(17),STAT) C C or new output table C ELSE IF ((APPFLG(1:1).EQ.'A') .OR. (APPFLG(2:2).EQ.'A')) THEN OUTFL = 4 CALL TBTOPN(OUTTAB,F_IO_MODE,TIDB,STAT) DO 150 IAV=1,OUTCOL !test, if label(s) already there CALL TBLSER(TIDB,LABEL(IAV),OUCOLN(IAV),STAT) IF (OUCOLN(IAV).LE.0) !Not found, error.. + CALL STETER(9,'Missing columns in output table...') 150 CONTINUE CALL TBLSER(TIDB,LABEL(17),OUCOLN(17),STAT) CALL TBIGET(TIDB,N,NAPP,N,N,N,STAT) !get rows already there C ELSE OUTFL = 3 !for new table, create it CALL TBTINI(OUTTAB,0,F_O_MODE,25,100,TIDB,STAT) DO 200 IAV=1,6 CALL TBCINI(TIDB,D_R8_FORMAT,1,'E24.15', + TABUNI(IAV),LABEL(IAV),OUCOLN(IAV),STAT) 200 CONTINUE DO 220 IAV=7,OUTCOL CALL TBCINI(TIDB,D_R4_FORMAT,1,'E12.5', + TABUNI(IAV),LABEL(IAV),OUCOLN(IAV),STAT) 220 CONTINUE C CALL TBCINI(TIDB,D_C_FORMAT,8,'A8', + TABUNI(17),LABEL(17),OUCOLN(17),STAT) ENDIF ENDIF C C +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C here all output options join: C display header line + branch according to input option C +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C 500 IF (METHOD(1:1).EQ.'I') THEN CBUF(1:) = ' xcenter ycenter angle'// + ' angl-err' CALL STTPUT(CBUF,STAT) CBUF(1:) = ' ' CBUF(35:) = 'maj-axis min-axis'// + ' icent ident' ELSE CBUF(1:) = ' xcenter ycenter xerr'// + ' yerr' CALL STTPUT(CBUF,STAT) CBUF(1:) = ' ' CBUF(35:) = ' xfwhm yfwhm'// + ' icent ident' ENDIF CALL STTPUT(CBUF,STAT) C C branch according to input option NMAL = 0 !init counter IF (INFLAG.EQ.1) THEN GOTO 900 ELSE IF (INFLAG.EQ.2) THEN GOTO 2000 ELSE IF (INFLAG.EQ.4) THEN GOTO 1500 ENDIF C C -------------------------------------------------------------------- C C C a) complete input frame is used - extract max. 200 * 200 window C around crude center C INFLAG = 3 (no looping) C C -------------------------------------------------------------------- C N = NPIX(1)/2 - 100 !get start of window IF (N.LE.0) N = 1 IF (ONEDIM.EQ.1) THEN NN = 1 ELSE NN = NPIX(2)/2 - 100 IF (NN.LE.0) NN = 1 ENDIF SPIX(1) = N SPIX(2) = NN C N = N + 199 !get end of window IF (N.GT.NPIX(1)) N = NPIX(1) IF (ONEDIM.EQ.1) THEN NN = 1 ELSE NN = NN + 199 IF (NN.GT.NPIX(2)) NN = NPIX(2) ENDIF SPIX(3) = N SPIX(4) = NN C DIN(1) = SPIX(1) DIN(2) = SPIX(2) CALL FPXWCO(1,IMNO,DIN,DOUT,STAT) !move from fp's to wco's DVALS(1) = DOUT(1) DVALS(2) = DOUT(2) DIN(1) = SPIX(3) DIN(2) = SPIX(4) CALL FPXWCO(1,IMNO,DIN,DOUT,STAT) DVALS(3) = DOUT(1) DVALS(4) = DOUT(2) DO 800,IAV=1,4 RVALS(IAV) = REAL(DVALS(IAV)) 800 CONTINUE NMAL = 1 !maybe we write a row to table GOTO 3000 C C -------------------------------------------------------------------- C C b) cursor defined subimage(s) - get main control block of ImageDisplay C INFLAG = 1 (loop on cursor input) C C -------------------------------------------------------------------- C 900 CALL DTOPEN(1,STAT) CALL STKRDI('DAZIN',1,1,IAV,INBUF,UNI,NULO,STAT) CALL STKRDI('DAZHOLD',14,1,IAV,OVCON,UNI,NULO,STAT) CALL STKRDI('CURSOR',1,3,IAV,INPUT,UNI,NULO,STAT) IF (INPUT(1).EQ.1) THEN NOCURS = 0 CALL STKRDI('INPUTI',1,2,IAV,NW,UNI,NULO,STAT) IF (NW(2).GT.MAPROW) THEN !free allocated space CALL STFYMP(-1,D_R4_FORMAT,TIMNO,PNTRA,STAT) MAPROW = NW(2) N = NPIX(1)*MAPROW CALL STFYMP(N,D_R4_FORMAT,TIMNO,PNTRA,STAT) ENDIF XU = NW(1) / 2 XL = XU IF ((XU+XL).EQ.NW(1)) XL = XL - 1 YU = NW(2) / 2 YL = YU IF ((YU+YL).EQ.NW(2)) YL = YL - 1 ELSE NOCURS = 2 ENDIF DRAWY(2) = INPUT(2) NLIM = INPUT(3) IF (NLIM.LT.1) NLIM = 1 C CALL DTGICH(QDSPNO,QIMCH,CBUF,RINF,STAT) CALL STKRDI('AUX_MODE',9,1,IAV,N,UNI,NULO,STAT) IF (N.NE.0) CALL CONCHA(QDSPNO,QOVCH,1,0) !clear overlay channel IF (OVCON.EQ.QIMCH) + CALL DAZZSC(QDSPNO,QOVCH,ZOOMX,SCROLX,SCROLY,STAT) FRAME(1:) = ' ' !that's needed for cursor interface C IF (INBUF(1).GT.0) THEN !we use the zoom window... IF (NOCURS.EQ.0) THEN CACT(5:7) = 'ZC0' !prepare zoom_window stuff for GETCUR ELSE CACT(5:7) = 'ZC2' ENDIF DRAWY(3) = 1 INBUF(1) = (QDSZX * 3) / 4 !use 3/4 of main display INBUF(2) = (QDSZY * 3) / 4 INBUF(3) = 0 !start at lower left of screen INBUF(4) = 0 CALL STKWRI('DAZIN',INBUF,2,4,UNI,STAT) N = 0 CALL SETCUR(QDSPNO,N,3,2,COOS,STAT) !use open_cross C C if `novice', display initial help, show current zoomfactor CALL STKRDI('ERROR',2,1,IAV,N,UNI,NULO,STAT) IF (N.EQ.1) CALL AUXHLP(0) !for novice users display help ELSE IF (NOCURS.EQ.0) THEN CALL SETCUR(QDSPNO,NOCURS,3,2,COOS,STAT) !open_cross ELSE CALL SETCUR(QDSPNO,NOCURS,1,2,COOS,STAT) !rectangle ENDIF ENDIF CALL PIXXCV('INIT',IMNO,RBUF,N) C C input cursor positions and check status 1000 IF (IDFLAG.EQ.0) THEN WRITE(IDF(3:),10010) NMAL+1 ELSE CALL STTPUT + ('Enter identifier (cursor must be in display window) :', + STAT) !ask for identifier... IAV = 8 CALL GETSTR(IDF,IAV) IF (IDF.EQ.' ') IDF = OLDIDF OLDIDF = IDF !save last identifier ENDIF GOTO 1100 C C invalid/bad cursor coords 1010 CALL STTDIS('Invalid cursor coords. ...',99,STAT) C C here the real thing... C 1100 CALL GETCUR(CACT,FRAME, + ICUR1,PCUR1(1),RVALS(1),RV,STAT1, + ICUR2,PCUR2(1),RVALS(3),RV,STAT2) IF (NOCURS.EQ.0) STAT2 = 0 C C check cursor status IF ( (STAT1.EQ.0) .AND. (STAT2.EQ.0) ) THEN IF (NMAL.EQ.0) THEN IF (COOFF.EQ.1) GOTO 9000 CALL STTPUT(INFO,STAT) COOFF = 1 FRAME(1:) = ' ' GOTO 1000 ELSE GOTO 9000 !cursor loop terminated ... ENDIF ENDIF C C show graph used for centering IF (NOCURS.GT.1) THEN C C do nothing on degenerated rectangles IF ((ICUR1(1).EQ.ICUR2(1)).AND.(ICUR1(2).EQ.ICUR2(2))) THEN CALL STTDIS('increase size of cursor rectangle ...',99,STAT) GOTO 1000 ENDIF SPIX(1) = NINT(PCUR1(1)) SPIX(2) = NINT(PCUR1(2)) SPIX(3) = NINT(PCUR2(1)) SPIX(4) = NINT(PCUR2(2)) BOX(1) = ICUR1(1) BOX(2) = ICUR1(2) BOX(3) = ICUR2(1) BOX(4) = ICUR2(2) ELSE NN = NINT(PCUR1(1)) SPIX(1) = NN - XL IF (SPIX(1).LT.1) SPIX(1) = 1 !check limits SPIX(3) = NN + XU IF (SPIX(3).GT.NPIX(1)) SPIX(3) = NPIX(1) NN = NINT(PCUR1(2)) SPIX(2) = NN - YL IF (SPIX(2).LT.1) SPIX(2) = 1 SPIX(4) = NN + YU IF (SPIX(4).GT.NPIX(2)) SPIX(4) = NPIX(2) PCUR1(1) = SPIX(1) !frame pixels of rectangle PCUR1(2) = SPIX(2) PCUR2(1) = SPIX(3) PCUR2(2) = SPIX(4) CALL PIXXCV('_RS',0,PCUR1,STAT) !convert to screen pixels IF (STAT.NE.0) GOTO 1010 CALL PIXXCV('_RS',0,PCUR2,STAT) IF (STAT.NE.0) GOTO 1010 BOX(1) = PCUR1(5) BOX(2) = PCUR1(6) BOX(3) = PCUR2(5) BOX(4) = PCUR2(6) CALL PIXXCV('_RW',0,PCUR1,STAT) !we also need world coordinates IF (STAT.NE.0) GOTO 1010 CALL PIXXCV('_RW',0,PCUR2,STAT) IF (STAT.NE.0) GOTO 1010 RVALS(1) = PCUR1(5) RVALS(2) = PCUR1(6) RVALS(3) = PCUR2(5) RVALS(4) = PCUR2(6) ENDIF C C get more precise wco's DIN(1) = PCUR1(1) DIN(2) = PCUR1(2) CALL FPXWCO(1,IMNO,DIN,DVALS(1),STAT) DIN(1) = PCUR2(1) DIN(2) = PCUR2(2) CALL FPXWCO(1,IMNO,DIN,DVALS(3),STAT) C C draw surrounding rectangle CALL DRAWCN(DRAWY,BOX,OLDBOX,STAT) IF (STAT.NE.0) THEN WRITE(CBUF,20002) NCOU CALL STTDIS(CBUF,99,STAT) NCOU = NCOU + 1 GOTO 1000 ENDIF C C continue at common section IF (SPIX(2).EQ.SPIX(4)) ONEDIM = 1 NMAL = NMAL + 1 !update coords. counter GOTO 3000 C C -------------------------------------------------------------------- C C C c) subimage around xfpix,yfpix C C C -------------------------------------------------------------------- C 1500 CALL STKRDI('INPUTI',1,2,IAV,NW,UNI,NULO,STAT) N = FRPIX(1) - NW(1)/2 !get start of window IF (N.LE.0) N = 1 IF (ONEDIM.EQ.1) THEN NN = 1 ELSE NN = FRPIX(2) - NW(2)/2 IF (NN.LE.0) NN = 1 ENDIF SPIX(1) = N SPIX(2) = NN C N = N + NW(1) !get end of window IF (N.GT.NPIX(1)) N = NPIX(1) IF (ONEDIM.EQ.1) THEN NN = 1 ELSE NN = NN + NW(2) IF (NN.GT.NPIX(2)) NN = NPIX(2) ENDIF SPIX(3) = N SPIX(4) = NN C DIN(1) = SPIX(1) DIN(2) = SPIX(2) CALL FPXWCO(1,IMNO,DIN,DOUT,STAT) !move from fp's to wco's IF (STAT.NE.0) THEN WRITE(CBUF,10070) DIN(1),DIN(2) CALL STTPUT(CBUF,STAT) ENDIF DVALS(1) = DOUT(1) DVALS(2) = DOUT(2) DIN(1) = SPIX(3) DIN(2) = SPIX(4) CALL FPXWCO(1,IMNO,DIN,DOUT,STAT) IF (STAT.NE.0) THEN WRITE(CBUF,10070) DIN(1),DIN(2) CALL STTPUT(CBUF,STAT) ENDIF DVALS(3) = DOUT(1) DVALS(4) = DOUT(2) DO 1800,IAV=1,4 RVALS(IAV) = REAL(DVALS(IAV)) 1800 CONTINUE NMAL = 1 !maybe we write a row to table GOTO 3000 C C -------------------------------------------------------------------- C C C d) table defined subimage(s) C C C -------------------------------------------------------------------- C 2000 NOCURS = 2 CALL TBIGET(TIDA,N,NROW,N,N,N,STAT) !get total no. of rows IF (NROW .LT. 1) CALL STETER (8,' Empty input table...') C CALL TBLSER(TIDA,LABEL(1),INCOLN(1),STAT) !look for :XSTART IF (INCOLN(1).LE.0) THEN CALL TBLSER(TIDA,GLABEL,INCOLN(1),STAT) !look for :X_POSITION IF (INCOLN(1).LE.0) THEN CBUF(1:) = 'Neither column :XSTART nor :X_POSITION '// + 'in input table found...' CALL STETER(4,CBUF) ELSE ONEDIM = 1 NOCURS = 1 !only 1 value (start/end) per row NROW = 2 * (NROW/2) !so make sure NROW is a multiple of 2 ENDIF C ELSE IF (ONEDIM.NE.1) THEN CALL TBLSER(TIDA,LABEL(2),INCOLN(2),STAT) !look for :YSTART ELSE CALL TBLSER(TIDA,LABEL(3),INCOLN(2),STAT) !look for :XEND ENDIF IF (INCOLN(2).LE.0) THEN CBUF(1:) = 'column labelled '//LABEL(IAV)// + ' is missing in input table...' CALL STETER(4,CBUF) ELSE IF (ONEDIM.NE.1) THEN DO 2040 IAV=3,4 !search for end columns CALL TBLSER(TIDA,LABEL(IAV),INCOLN(IAV),STAT) IF (INCOLN(IAV).LE.0) THEN CBUF(1:) = 'column labelled '//LABEL(IAV)// + ' is missing in input table...' CALL STETER(4,CBUF) ENDIF 2040 CONTINUE ENDIF ENDIF C C look for IDENT column 2080 CALL TBLSER(TIDA,LABEL(17),INCOLN(17),STAT) IF (INCOLN(17).NE.-1) THEN IDYES = 1 ELSE IDYES = 0 IF (OUTFL.EQ.2) CALL TBCINI(TIDA,D_C_FORMAT,8,'A8', + TABUNI(17),LABEL(17),OUCOLN(17),STAT) ENDIF C C here we loop 2100 NMAL = NMAL + 1 IF (NMAL.GT.NROW) THEN !test for end of table WRITE(CBUF,10001) NCOU,' table entries processed ' CALL STTPUT(CBUF,STAT) GOTO 9000 ENDIF C C get next row of values + read as double data CALL TBSGET(TIDA,NMAL,SELFLG,STAT) IF (.NOT.SELFLG) GOTO 2100 C IF (IDYES.EQ.1) THEN CALL TBERDC(TIDA,NMAL,INCOLN(17),IDF,TABNUL,STAT) ELSE WRITE(IDF(3:),10010) NMAL ENDIF C IF (ONEDIM.NE.1) THEN CALL TBRRDD(TIDA,NMAL,4,INCOLN,DVALS(1),TABNUL,STAT) DIN(1) = DVALS(1) !convert from wco's to fp's DIN(2) = DVALS(2) CALL FPXWCO(-1,IMNO,DIN,DOUT,STAT) IF (STAT.NE.0) THEN WRITE(CBUF,10077) NMAL CALL STTPUT(CBUF,STAT) GOTO 2100 ENDIF SPIX(1) = NINT(DOUT(1)) SPIX(2) = NINT(DOUT(2)) DIN(1) = DVALS(3) DIN(2) = DVALS(4) CALL FPXWCO(-1,IMNO,DIN,DOUT,STAT) IF (STAT.NE.0) THEN WRITE(CBUF,10077) NMAL CALL STTPUT(CBUF,STAT) GOTO 2100 ENDIF SPIX(3) = NINT(DOUT(1)) SPIX(4) = NINT(DOUT(2)) ELSE IF (NOCURS.EQ.1) THEN CALL TBRRDD(TIDA,NMAL,1,INCOLN,DVALS(1),TABNUL,STAT) DIN(1) = DVALS(1) CALL FPXWCO(-1,IMNO,DIN,DOUT,STAT) IF (STAT.NE.0) THEN WRITE(CBUF,10077) NMAL CALL STTPUT(CBUF,STAT) GOTO 2100 ENDIF SPIX(1) = NINT(DOUT(1)) NMAL = NMAL + 1 CALL TBRRDD(TIDA,NMAL,1,INCOLN,DVALS(3),TABNUL,STAT) DIN(1) = DVALS(3) CALL FPXWCO(-1,IMNO,DIN,DOUT,STAT) IF (STAT.NE.0) THEN WRITE(CBUF,10077) NMAL CALL STTPUT(CBUF,STAT) GOTO 2100 ENDIF SPIX(3) = NINT(DOUT(1)) ELSE CALL TBRRDD(TIDA,NMAL,2,INCOLN,DVALS(1),TABNUL,STAT) DVALS(3) = DVALS(2) DIN(1) = DVALS(1) CALL FPXWCO(-1,IMNO,DIN,DOUT,STAT) IF (STAT.NE.0) THEN WRITE(CBUF,10077) NMAL CALL STTPUT(CBUF,STAT) GOTO 2100 ENDIF SPIX(1) = NINT(DOUT(1)) DIN(1) = DVALS(3) CALL FPXWCO(-1,IMNO,DIN,DOUT,STAT) IF (STAT.NE.0) THEN WRITE(CBUF,10077) NMAL CALL STTPUT(CBUF,STAT) GOTO 2100 ENDIF SPIX(3) = NINT(DOUT(1)) ENDIF DVALS(2) = 0.0 !put into 2-dimj order DVALS(4) = 0.0 SPIX(2) = 1 SPIX(4) = 1 ENDIF C DO 2200,IAV=1,4 !save also in real array RVALS(IAV) = REAL(DVALS(IAV)) 2200 CONTINUE NCOU = NCOU + 1 C C +++++++++++++++++++++++++++++++++++++ C C common section for option a) and b) and c) C C +++++++++++++++++++++++++++++++++++++ C C check subimages dimensions 3000 NPIXT(1) = SPIX(3) - SPIX(1) + 1 !get size of subimage NPIXT(2) = SPIX(4) - SPIX(2) + 1 XOFF = SPIX(1) - 1. YOFF = SPIX(2) - 1. INBUF(1) = NPIX(1) INBUF(2) = NPIXT(2) INBUF(3) = 1 C IF (NPIXT(1).LT.3) THEN CBUF(1:) = 'x-dimension < 3 pixels -'// + ' subimage skipped...' CALL STTPUT(CBUF,STAT) GOTO 7000 ELSE IF (NPIXT(2) .GT. MAPROW) THEN CALL STFYMP(-1,D_R4_FORMAT,TIMNO,PNTRA,STAT) MAPROW = NPIXT(2) N = NPIX(1)*MAPROW CALL STFYMP(N,D_R4_FORMAT,TIMNO,PNTRA,STAT) ENDIF IF ((NPIXT(1)*NPIXT(2)).GT.40000) THEN !200 * 200 max. subwindow CBUF(1:) = 'subimage > 200*200 pixels - subimage skipped... ' CALL STTPUT(CBUF,STAT) GOTO 7000 ENDIF C C get a y-stripe from disk + copy subwindow out N = INBUF(1) * INBUF(2) !size of y-dim stripe in frame NN = ((SPIX(2)-1)*NPIX(1)) + 1 !1. pixel in there IF ((ONEDIM.NE.1).OR.(FIRTIM.EQ.0)) + CALL STFGET(IMNO,NN,N,IAV,MADRID(PNTRA),STAT) FIRTIM = 1 C ALOW(1) = SPIX(1) !first x,y pixel in source ALOW(2) = 1 !and destination frame HIGH(1) = SPIX(3) HIGH(2) = NPIXT(2) CALL COPWND(MADRID(PNTRA),INBUF(1),MADRID(PNTRC),NPIXT, + ALOW,BLOW,HIGH) IF (INVFLG.EQ.1) + CALL UPSIDE(MADRID(PNTRC),NPIXT,XCUTS) !move to `normal' Gaussian C C actual centering algorithm C CPIX(1) = 1 !xstart CPIX(2) = NPIXT(1) !xend CPIX(3) = 1 !ystart CPIX(4) = NPIXT(2) !yend C IF (METHOD(1:1) .EQ. 'I') THEN CALL IQEFUN(MADRID(PNTRC),NPIXT,RBUF(1),RBUF(11),KSTAT) RV = 1.0 !added to offset later on XOUT = RBUF(1) XSIG = RBUF(2)/2.35482 XERR = RBUF(11) YOUT = RBUF(3) YSIG = RBUF(4)/2.35482 YERR = RBUF(13) ICENT = RBUF(6) ANGLE = RBUF(5) ANGSIG = RBUF(15) IF (XSIG .GE. YSIG) THEN !the larger one -> major axis AXMAJ = RBUF(2)*STEP(1) AXMIN = RBUF(4)*STEP(2) ELSE AXMAJ = RBUF(4)*STEP(2) AXMIN = RBUF(2)*STEP(1) ENDIF ELSE CALL STACEN(MADRID(PNTRC),NPIXT(1),NPIXT(2),METHOD,CPIX, + XOUT,YOUT,XERR,YERR,XSIG,YSIG,ICENT,KSTAT) RV = 0.0 ANGLE = 0.0 ANGSIG = 0.0 ENDIF NCOOS = NCOOS + 1 !increment coord. counter C C check status XOUT = XOUT + XOFF + RV !relate to complete frame IF (KSTAT.EQ.0) THEN CSTAT(1:) = ' ' XERR = XERR*ABS( STEP(1) ) XSIG = XSIG*ABS( STEP(1) ) XFWHM = XSIG*2.35482 ELSE XERR = 0. XSIG = 0. XFWHM = 0. IF (KSTAT.LT.0) THEN !must be from IQEFUN CSTAT(1:) = 'estimation failed ' ELSE IF (KSTAT.EQ.1) THEN CSTAT(1:) = 'ambiguos sources ' ELSE IF (KSTAT.EQ.2) THEN CSTAT(1:) = 'no source ' ELSE IF (KSTAT.EQ.3) THEN CSTAT(1:) = 'iteration failed ' ELSE CSTAT(1:) = 'routine error ' ENDIF ENDIF C DIN(1) = XOUT IF (NPIXT(2).GT.1) THEN YOUT = YOUT + YOFF + RV !relate to complete frame IF (KSTAT.EQ.0) THEN YERR = YERR*ABS( STEP(2) ) YSIG = YSIG*ABS( STEP(2) ) YFWHM = YSIG*2.35482 ELSE YERR = 0. YSIG = 0. YFWHM = 0. ENDIF DIN(2) = YOUT ELSE YERR = 0. YSIG = 0. YFWHM = 0. DIN(2) = SPIX(2) ENDIF CALL FPXWCO(1,IMNO,DIN,DVALS(5),STAT) XOUT = REAL(DVALS(5)) !DVALS(1 - 4) hold are coords YOUT = REAL(DVALS(6)) C C for cursor input draw numbers... IF ( (INFLAG.EQ.1) .AND. (DRAWY(2).GT.1) ) THEN IF (IDFLAG.EQ.1) THEN LABL = IDF NN = 1 LLABL = INDEX(LABL,' ') - 1 IF (LLABL.LE.0) LLABL = LEN(LABL) ELSE WRITE(LABL(5:),10011) NCOOS NN = 5 DO 3300 N=5,7 IF (LABL(N:N).EQ.' ') NN = N + 1 3300 CONTINUE LLABL = 9 - NN ENDIF C XYALP(2) = ICUR1(2) - 9 !place string below N = ICUR2(1) - ICUR1(1) LLABL = 7 * LLABL IF (N.LE.LLABL) THEN XYALP(1) = ICUR1(1) ELSE N = (N-LLABL) / 2 XYALP(1) = ICUR1(1) + N ENDIF CALL IIGTXT(QDSPNO,QOVCH,LABL(NN:),XYALP(1),XYALP(2), + 0,0,255,1,STAT) ENDIF C C Before we display the results, we should check if FWHM values make sense IF ( ( XFWHM .GT. (NPIXT(1) * ABS(STEP(1)) ) ) .OR. + ( YFWHM .GT. (NPIXT(2) * ABS(STEP(2)) ) ) ) THEN KSTAT = 4 CSTAT(1:) = 'FWHM exceeds size of image ' ENDIF IF ((XFWHM .LT. 0.0) .OR. (YFWHM .LT. 0.0)) THEN KSTAT = 5 CSTAT(1:) = 'Negative FWHM calculated ' ENDIF C C display results RVALS(5) = XOUT RVALS(6) = YOUT RVALS(7) = XERR RVALS(8) = YERR RVALS(9) = XSIG RVALS(10) = YSIG RVALS(11) = KSTAT RVALS(14) = ICENT C IAV = 14 IF (CSTAT(1:1).EQ.' ') THEN IF (METHOD(1:1) .EQ.'I') THEN RVALS(12) = AXMAJ RVALS(13) = AXMIN RVALS(15) = ANGLE RVALS(16) = ANGSIG WRITE(CBUF,20100) RVALS(5),RVALS(6),RVALS(15),RVALS(16) CALL STTPUT(CBUF,STAT) WRITE(CBUF,20101) RVALS(12),RVALS(13),RVALS(14),IDF IAV = 16 ELSE RVALS(12) = XFWHM RVALS(13) = YFWHM WRITE(CBUF,20100) RVALS(5),RVALS(6),RVALS(7),RVALS(8) CALL STTPUT(CBUF,STAT) WRITE(CBUF,20101) RVALS(12),RVALS(13),RVALS(14),IDF ENDIF ELSE RVALS(12) = 0. RVALS(13) = 0. WRITE(CBUF,20020) RVALS(5),RVALS(6),CSTAT,IDF ENDIF CALL STTPUT(CBUF,STAT) CBUF = ' ' !clear last line CALL STTPUT(CBUF,STAT) C CALL STKWRR('OUTPUTR',RVALS(1),1,IAV,UNI,STAT) CALL STKWRD('OUTPUTD',DVALS(1),1,6,UNI,STAT) C C fill descriptor, table or nothing IF (OUTFL.EQ.1) THEN CALL STDWRR(IMNO,DESCR,RVALS,FELEM,IAV,UNI,STAT) !fill descr FELEM = -1 C ELSE IF (OUTFL.GT.1) THEN IF (OUTFL.EQ.2) THEN !same table IAV = 5 NN = NMAL !use same rows ELSE !different table IAV = 1 IF (OUTFL.EQ.4) THEN NAPP = NAPP + 1 NN = NAPP !append new rows ELSE !use sequential, new rows IF (INFLAG.NE.2) THEN !no table input NN = NMAL ELSE NN = NCOU !table input ENDIF ENDIF ENDIF YL = 7 - IAV !write 6 or 2 double columns CALL TBRWRD(TIDB,NN,YL,OUCOLN(IAV),DVALS(IAV),STAT) IAV = 7 CALL TBRWRR(TIDB,NN,4,OUCOLN(IAV),RVALS(IAV),STAT) RV = KSTAT CALL TBRWRR(TIDB,NN,1,OUCOLN(14),RV,STAT) CALL TBRWRR(TIDB,NN,3,OUCOLN(11),RVALS(12),STAT) IF (METHOD(1:1) .EQ.'I') + CALL TBRWRR(TIDB,NN,2,OUCOLN(15),RVALS(15),STAT) CALL TBEWRC(TIDB,NN,OUCOLN(17),IDF,STAT) ENDIF C C and loop again 7000 IF (INFLAG.EQ.1) THEN IF (NMAL.LT.NLIM) GOTO 1000 ELSE IF (INFLAG.EQ.2) THEN GOTO 2100 ENDIF C C That's it folks... 9000 IF (OUTFL.EQ.3) CALL TBSINI(TIDB,STAT) IF (TIDB.GE.0) THEN IF (ACTION(1:1).EQ.'G') THEN IDENT(1:) = 'GAUSS ' N = 5 ELSE IF (ACTION(1:1).EQ.'I') THEN IDENT(1:) = 'IQE ' N = 3 ELSE IDENT(1:) = 'MOMENT ' N = 6 ENDIF CBUF(1:) = + 'From CENTER/'//IDENT(1:N)//' using frame: '//FRAME CALL STDWRC(TIDB,'HISTORY',1,CBUF,-1,80,UNI,STAT) ENDIF C IF (INFLAG.EQ.1) THEN CALL DTCLOS(QDSPNO) CALL REFOVR(STAT) !refresh overlay ENDIF C C save no. of coordinates obtained for subsequent applications CALL STKWRI('OUTPUTI',NCOOS,1,1,UNI,STAT) CALL STSEPI C C format statements 10001 FORMAT(I4,A) 10010 FORMAT(I4.4) 10011 FORMAT(I4) 10070 FORMAT('bad coord(s) - X, Y: ',2F20.10) 10077 FORMAT('row no. ',I5,' contains bad coord(s) - we skip ...') 20000 FORMAT(G10.4,1X,G10.4,1X,G8.2,1X,G8.2,1X,G9.3,1X,G9.3, + 1X,G10.4,1X,A8) 20002 FORMAT(I4,' cursor rectangle too small...') 20020 FORMAT(G10.4,1X,G10.4,10X,A,8X,A) 20100 FORMAT(G14.8,1X,G14.8,1X,G12.5,2X,G12.5) 20101 FORMAT(32X,G12.5,1X,G12.5,2X,G12.5,1X,A8) C END SUBROUTINE UPSIDE(A,MPIX,RMAX) C IMPLICIT NONE C REAL A(*),RMAX C INTEGER MPIX(2) INTEGER N, OFF C OFF = 1 DO 1000,N=1,MPIX(1)*MPIX(2) A(OFF) = RMAX - A(OFF) OFF = OFF + 1 1000 CONTINUE C RETURN END SUBROUTINE DRAWCN(DRAWY,NBOX,OLDBOX,STAT) C IMPLICIT NONE C INTEGER DRAWY(3),NBOX(4),OLDBOX(4),STAT INTEGER INFO(11) INTEGER BOX(4) INTEGER OFIGX(5),OFIGY(5) INTEGER OFAGX(5),OFAGY(5) INTEGER NPO,NAPO INTEGER DSPNO,OVL,ZOOM,N INTEGER XOFF,YOFF INTEGER SIZE,SIZEA C REAL RTEMP(2) C INCLUDE 'MID_INCLUDE:IDIDEV.INC' INCLUDE 'MID_INCLUDE:IDIMEM.INC' C DATA RTEMP /-1.0,-1.0/ C C test DRAWY options STAT = 0 C C here for start up - update disp no. + overlay, if we use aux_window IF (DRAWY(1).EQ.-1) THEN DRAWY(1) = 0 IF (DRAWY(3).EQ.1) THEN CALL AUXWND(6,INFO,BOX,BOX,STAT) DSPNO = INFO(1) !get display no. of aux_window OVL = 1 ELSE DSPNO = QDSPNO !use main window OVL = QOVCH ZOOM = 1 XOFF = 0 YOFF = 0 ENDIF C C here for DRAWY(1) = 0 ELSE DO 200 N=1,4 !check, if cursor moved IF (NBOX(N).NE.OLDBOX(N)) GOTO 1000 200 CONTINUE RETURN !nothing to do... ENDIF C C test, if we use aux_window 1000 IF (DRAWY(3).EQ.1) THEN CALL IIMCMY(DSPNO,OVL,1,0,STAT) CALL AUXWND(6,INFO,BOX,BOX,STAT) ZOOM = INFO(4) XOFF = INFO(6) !lower left screen pixel of YOFF = INFO(7) !extracted square ENDIF C C get center and size of new box SIZE = NBOX(3) - NBOX(1) + 1 SIZEA = NBOX(4) - NBOX(2) + 1 IF (SIZE.GT.SIZEA) SIZE = SIZEA !take min of x-, y-size C C build + draw box for aux_window + main display IF (DRAWY(3).EQ.1) THEN BOX(1) = (NBOX(1) - XOFF) * ZOOM !take out offset BOX(2) = (NBOX(2) - YOFF) * ZOOM !of extracted zoom square BOX(3) = (NBOX(3) - XOFF) * ZOOM !and BOX(4) = (NBOX(4) - YOFF) * ZOOM !adapt to zoom factor CALL BLDGRA('REC',BOX,RTEMP,OFAGX,OFAGY,5,NAPO) CALL IIGPLY(DSPNO,OVL,OFAGX,OFAGY,NAPO,255,1,STAT) ENDIF IF (DRAWY(2).EQ.1) THEN !box for original channel CALL BLDGRA('REC',NBOX,RTEMP,OFIGX,OFIGY,5,NPO) CALL IIGPLY(QDSPNO,QOVCH,OFIGX,OFIGY,NPO,255,1,STAT) ENDIF C C finally update OLDBOX DO 5000, N=1,4 OLDBOX(N) = NBOX(N) 5000 CONTINUE C C that's it folks... C STAT = 0 RETURN END