C @(#)pos1a.for 17.1.1.1 (ESO-DMD) 01/25/02 17:13:38 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 POS1a *--------------------------------------------------------------------------- c.file: pos1a.for c.purpose: POS1 astrometric package C THIS PROGRAM COMPUTES (ALFA,DELTA) FOR OBJECTS C MEASURED ON PHOTOGRAPHIC PLATES. C C.remarks: C The plates can be from the Schmidt or any other telescope. If they c are from the Schmidt, the plate curvature will be taken into account. c c INPUT is from two midas tables: c STANDARD: the following columns must exist c :PPM Identifier (<=6digit intg) c :R_A Right Ascension in deg. c :DEC Decl. in deg c :MAG Magnitude c :PMA Proper Motion in RA in arcsecond/year (NOT time_second) c :PMD Proper Motion in Dec in arcsecond/year C DATA FILE: the following columns must exist c :IDENT Identifier; must contain a number, leading c letters (eg, PPM) are ignored c :Xcen X measurement c :YCEN Y Measurement C C The following parameters are read from MIDAS keywords: C FLAGS(1): Schmidt plate? (takes the curvature into account) C FLAGS(2): Measured with the La Silla blink? c SIG: Permitted measurement RMS, in measurement unit (eg micron) C CENTCOOR: coordinates of the plate center; either C h,m,s,d,'," (the routine take care of the -0deg), or C MEAN: takes the barycenter of standard stars C EPO_PLA: epoch of the plate C EPO_CAT: epoch of the catalogue (NB: This is NOT C necessaryly the equinox of the catalogue coordinates. C TERMS: which from the following terms should be used C for the transfo. 1: set, 0: ignored; same for Y C X Y XY X**2 Y**2 X**3 Y**3 XY**2 YX**2 C STDSEL: Standard selection C A: take ALL the available standards (eg, 1st run) C U: take only the standard not marked as deleted (eg, C subsequent run). C OUTPUT: C The parameter of the transformation are stored in the C Following Midas keywords: C CX, CY: constant terms C BX, BY: parameters of x and y and their power C NX, NY: how many terms in x and y c ALDE0: Alpha, Delta, cos Delta, sin Delta C C Call this prg with the pos1.prg or pos1a.prg. C pos1a.f computes ONLY the transformation parameters. C The actual coordinate conversion is performed by pos1b. C The standard can be selected/deleted by residual.prg *--------------------------------------------------------------------------- c.author: OR Hainaut C.versions: C - GENEVA, DECEMBER 16, 1976 R.M.West C - REVISED SEPTEMBER 20, 1977 RMW C - MUNICH, THOROUGHLY REVISED OCTOBER 1980 RMW C - TRANSFERED TO VAX 31.1.1985 RMW C - INSTALLED AT LA SILLA 22.7.1985 RMW C - REVISED MARCH 1987 (PROPER MOTION R.A.) RMW C - transfered to Sun/Unix at la Silla, 1992, A Smette & OR Hainaut C - re-installed on Sun/Unix at Garching, 9/92 ORH C - Structure cleaned (IF-THEN-ELSE, GOTO...) 5/93 ORH c - transfered to MIDAS 8/93 ORH c - a bug fixed: accuracy similar to original POS1 C #Version: Sat Sep 17 19:18:47 1994 C.VERSION 1995-Mar-26 : Use PARAMETER + Increase buffer size to 10000, PJG C C 000609 *---------------------------------------------------------------------------- C IMPLICIT REAL*8 (A-H,O-Z), INTEGER*4 (I-N) C PARAMETER (MAXOBJ=100001) C add parameter MAXSTD to increase the size PN12/00: PARAMETER (MAXSTD=100001) C CHARACTER*1 STDSEL,DSIG CHARACTER*2 ITEGN1,ITEGN2,FLAGS,KTEXT(3,9) CHARACTER*16 PSFILE,DFILE CHARACTER*16 ACUR, AFOR CHARACTER*19 STRING_TERM CHARACTER*72 CENTCOOR CHARACTER*250 MSG C here the variables for measurements INTEGER TIDD INTEGER NCOLD INTEGER NROWD INTEGER NCDIDENT INTEGER NSC,ACOL,AROW,ISTAT INTEGER ILOOP C here the variables for standards file: INTEGER TIDS INTEGER NCOLS INTEGER NROWS INTEGER NCSPPM LOGICAL LOG, NULL, SELECT, ROWSEL, SCHMIDT,BLINK, CENTER DIMENSION MNOST(MAXOBJ),IFOUND(MAXOBJ),NOST(MAXOBJ) DIMENSION KK(9),X(10),AR(MAXOBJ,4),BX(9),BY(9),IDX(9),IDY(9) DIMENSION A0(3),D0(3),ST(MAXSTD,2),XY(MAXOBJ,2),ALFA(3),DELTA(3) DIMENSION SSX(10,10),S(10,10),R(10,10),SIY(10),SLSSC(9) DIMENSION SDC(9),PRXY(MAXOBJ,2) DIMENSION VBI(9),B(9),SV(10),SIGMA(10),XBAR(10),SP(10) DIMENSION SSINV(9,9),ERXY(3),PRST(MAXSTD,2) DIMENSION ISTDR(MAXOBJ), ISTDR2(MAXSTD) DIMENSION XTEMP(20),YTEMP(20),RESA(3),RESB(3) DIMENSION STMAG(MAXSTD),IA1(2),ID1(2),IA2(2),ID2(2) DIMENSION KICK(MAXSTD),PRSTMAG(MAXSTD) EQUIVALENCE > (SSX,S,R), > (S(2,1),B), > (S(3,2),DEV), > (S(4,2),SW,V), > (S(5,2),SSR,RSW), > (S(6,2),CONST), > (S(7,2),SSDR,SUM), > (S(8,2),RSQD), > (S(9,2),CHG,ROFF), > (S(4,3),XMSDR), > (S(5,3),XMSDD), > (S(6,3),VE), > (S(7,3),SDE), > (S(1,4),SLSSC), > (S(1,5),TSQDC), > (S(1,6),SDC), > (S(1,7),VBI), > (SV,SIGMA), > (X,XBAR) *--------------------------------------------------------------------- *--- MIDAS... *--------------------------------------------------------------------- INCLUDE 'MID_INCLUDE:ST_DEF.INC' COMMON /VMR/MADRID INCLUDE 'MID_INCLUDE:ST_DAT.INC' DATA KTEXT/' X', ' ', ' ', ' Y', ' ', ' ', ' X', 'Y ', ' ', > ' X', '**', '2 ', ' Y', '**', '2 ', ' X', '**', '3 ', > ' Y', '**', '3 ', ' X', 'Y*', '*2', ' Y', 'X*', '*2'/ *--------------------------------------------------------------------- *--- init CALL STSPRO('POS1A') *--------------------------------------------------------------------- LU = 30 DO0 = 0.0D0 DO1 = 1.0D0 DO12 = 1.2D1 DO24 = 2.4D1 DO60 = 6.0D1 DO180 = 1.8D2 DO3600 = 3.6D3 PI = 4.0D0*DATAN(1.0D0) DOPI18 = PI/DO180 FOC = 3.054D2/3.6D3 C PSFILE = ' ' DFILE = ' ' C DO 5, J=1,100 5 KICK(J) = 0 NKICK = 0 ITRY = 0 10 WRITE(MSG,9015) 9015 FORMAT(' POS1: Astrometric solution for Schmid/CCD images') CALL STTPUT(MSG,ISTAT) CALL STTPUT(' ',ISTAT) *--- Read the instrument flags: CALL STKRDC('FLAGS',1,1,2,N1,FLAGS,IDUM1,IDUM2,ISTAT) IF (FLAGS(1:1) .EQ. 'S' .OR. FLAGS(1:1) .EQ. 's' $ .OR. FLAGS(1:1).EQ.'Y' .OR. FLAGS(1:1).EQ.'y') THEN SCHMIDT = .TRUE. ELSE SCHMIDT = .FALSE. ENDIF IF (FLAGS(2:2) .EQ. 'y' .OR. FLAGS(2:2) .EQ. 'Y' > .OR. FLAGS(2:2) .EQ. 'B' .OR. FLAGS(2:2) .EQ. 'b' ) THEN BLINK = .TRUE. ELSE BLINK = .FALSE. ENDIF *--- Get the permitted permitted RMS CALL STKRDD('SIG',1,1,N1,SIG,IDUM1,IDUM2,ISTAT) *--- Get the center of plate CALL STKRDC('CENTCOOR',1,1,72,N1,CENTCOOR,IDUM1,IDUM2,ISTAT) IF (CENTCOOR.NE.'MEAN' .AND. CENTCOOR.NE.'mean') THEN * * Not "mean": try to read the coordinates. * * check if there is a '-': negative decl. RH = DO1 IF(INDEX(CENTCOOR,'-').GT.0) RH = -DO1 * * Reads the coord. READ(CENTCOOR,*,ERR=105,END=105) A0,D0 DO 103 I = 1, 3 103 D0(I) = DABS(D0(I)) ALFA0 = (A0(1)+A0(2)/DO60+A0(3)/DO3600)*1.5D1*DOPI18 DELTA0 = RH*(D0(1)+D0(2)/DO60+D0(3)/DO3600)*DOPI18 CENTER = .TRUE. GOTO 106 105 CENTER = .FALSE. 106 CONTINUE ELSE * * "Mean", or error reading the coordinates CENTER = .FALSE. ENDIF *--- Read the plate and catalogue epochs: CALL STKRDD('EPO_PLA',1,1,N1,EPOCH,IDUM1,IDUM2,ISTAT) CALL STKRDD('EPO_CAT',1,1,N1,EP1,IDUM1,IDUM2,ISTAT) *--- Reads the xy transformation terms selection flags CALL STKRDC('TERMS',1,1,19,N1,STRING_TERM,IDUM1,IDUM2,ISTAT) *--- read the standard stars status: * * A= check all standards * * U= check only the previously identified/undeleted stars. CALL STKRDC('STDSEL',1,1,1,N1,STDSEL,IDUM1,IDUM2,ISTAT) SELECT = .FALSE. IF (STDSEL .eq. 'U') SELECT = .TRUE. *--- initialize the ASCII log file cc LU_LOG = 30 cc OPEN(UNIT=LU_LOG,FILE='pos1.log') IF (SCHMIDT)THEN WRITE(MSG,9053) CALL STTPUT(MSG,ISTAT) ELSE WRITE(MSG,9153) CALL STTPUT(MSG,ISTAT) ENDIF 9053 FORMAT(' The plate is assumed to be a SCHMIDT plate (curved)') 9153 FORMAT(' The frame is assumed NOT to be a SCHMIDT', $ ' plate (NOT curved)') MSG = ' ' IF (BLINK) WRITE(MSG,9052) CALL STTPUT(MSG,ISTAT) 9052 FORMAT(' (X,Y) measurements in data file were measured by BLINK') *--------------------------------------------------------------------- *--- open the input files *--------------------------------------------------------------------- *--- *standards file CALL STKRDC('IN_STD',1,1,20,N1,PSFILE,IDUM1,IDUM2,ISTAT) * * opening the table CALL TBTOPN(PSFILE,F_IO_MODE,TIDS,ISTAT) * * check the numbers of rows and columns of the table CALL TBIGET(TIDS,NCOLS,NROWS,NSC,ACOL,AROW,ISTAT) * * check if the columns exists * * PPM identifier CALL TBCSER(TIDS,':PPM',NCSPPM,ISTAT) IF (NCSPPM.LE.0) THEN CALL STTPUT('*error* Input column PPM not found',ISTAT) STOP ENDIF * * magnitude CALL TBCSER(TIDS,':MAG',NCSMAG,ISTAT) IF (NCSMAG.LE.0) THEN CALL STTPUT('*error* Input column MAG not found',ISTAT) STOP ENDIF * * right asc. CALL TBCSER(TIDS,':R_A',NCSR_A,ISTAT) IF (NCSR_A.LE.0) THEN CALL STTPUT('*error* Input column R_A not found',ISTAT) STOP ENDIF * * declination CALL TBCSER(TIDS,':DEC',NCSDEC,ISTAT) IF (NCSDEC.LE.0) THEN CALL STTPUT('*error* Input column DEC not found',ISTAT) STOP ENDIF * * Proper motion in right asc. CALL TBCSER(TIDS,':PMA',NCSPMA,ISTAT) IF (NCSPMA.LE.0) THEN CALL STTPUT('*error* Input column PMA not found',ISTAT) STOP ENDIF * * proper motion in decl. CALL TBCSER(TIDS,':PMD',NCSPMD,ISTAT) IF (NCSPMD.LE.0) THEN CALL STTPUT('*error* Input column PMD not found',ISTAT) STOP ENDIF * * status CALL TBCSER(TIDS,':STD',NCSSTD,ISTAT) IF (NCSSTD.LE.0) THEN * * if the column does not exist, it's created. CALL STTPUT('*info* Input column :STD not found; created',ISTAT) CALL TBCINI( TIDS, D_I4_FORMAT, 1, 'I1', 'Flg', > 'STD', NCSSTD, ISTAT ) ENDIF * * error in x CALL TBCSER(TIDS,':XERR',NCSXER,ISTAT) IF (NCSXER.LE.0) THEN * * if the column does not exist, it's created. CALL STTPUT('*info* Input column XERR not found; created',ISTAT) CALL TBCINI( TIDS, D_R8_FORMAT, 1, 'F8.4' , > 'ArcSec', 'XERR', NCSXER ,ISTAT) ENDIF * * error in y CALL TBCSER(TIDS,':YERR',NCSYER,ISTAT) IF (NCSYER.LE.0) THEN * * if the column does not exist, it's created. CALL STTPUT('*info* Input column YERR not found; created',ISTAT) CALL TBCINI( TIDS, D_R8_FORMAT, 1, 'F8.4' , > 'ArcSec', 'YERR', NCSYER ,ISTAT) ENDIF *--- *measurements file CALL STKRDC('IN_MES',1,1,20,N1,DFILE,IDUM1,IDUM2,ISTAT) * * opening the table CALL TBTOPN(DFILE,F_I_MODE,TIDD,ISTAT) * * check the numbers of rows and columns of the table CALL TBIGET(TIDD,NCOLD,NROWD,NSC,ACOL,AROW,ISTAT) * * check if the columns exists * * PPM identifier CALL TBCSER(TIDD,':IDENT',NCDIDENT,ISTAT) IF (NCDIDENT.LE.0) THEN CALL STTPUT('*error* Input column IDENT not found',ISTAT) STOP ENDIF * * x position CALL TBCSER(TIDD,':XCEN',NCDX,ISTAT) IF (NCDX.LE.0) THEN CALL STTPUT('*error* Input column XCEN not found',ISTAT) STOP ENDIF * * y position CALL TBCSER(TIDD,':YCEN',NCDY,ISTAT) IF (NCDY.LE.0) THEN CALL STTPUT('*error* Input column YCEN not found',ISTAT) STOP ENDIF 55 DO I = 1,3 RESA(I) = 0.0 RESB(I) = 0.0 ENDDO ITOT = 0 IRES = 0 IRESNO = 0 *-- Write to the log file: WRITE(MSG,9056) PSFILE CALL STTPUT(MSG,ISTAT) 9056 FORMAT(' Standard stars from table : ',A16) WRITE(MSG,9156) DFILE CALL STTPUT(MSG,ISTAT) 9156 FORMAT(' XY measurements from table : ',A16) CALL STTPUT(' ',ISTAT) *--------------------------------------------------------------------- *--- Reads the measurements (X,Y) *--------------------------------------------------------------------- WRITE(MSG,9157) 9157 FORMAT(' --- (X,Y) MEASUREMENTS ---') CALL STTPUT(MSG,ISTAT) IF (BLINK) WRITE(MSG,9057) CALL STTPUT(MSG,ISTAT) 9057 FORMAT(' The following (X,Y)s are expressed in microns') WRITE(MSG,9058) SIG CALL STTPUT(MSG,ISTAT) 9058 FORMAT(' Mean measured (X,Y) after cleaning R.M.S. > ',F5.2) CALL STTPUT(' ',ISTAT) WRITE(MSG,9158) CALL STTPUT(MSG,ISTAT) 9158 FORMAT(' ID.',7X,'X',11X,'Y',7X,'N SigmaX SigmaY ', >'SigmaR Deleted') WRITE(MSG,9258) 9258 FORMAT(' <----> <---------> <---------> <->', $ ' <-----> <----> <----> <->') CALL STTPUT(MSG,ISTAT) C loop from ILOOP = 1 to ILOOP = number of measurement rows DO 62 ILOOP = 1, NROWD * *identifier CALL TBERDC(TIDD,ILOOP,NCDIDENT,ACUR,NULL,ISTAT) * *(x/ycen are r*4 read as r*8) CALL TBERDD(TIDD,ILOOP,NCDX,XCUR,NULL,ISTAT) CALL TBERDD(TIDD,ILOOP,NCDY,YCUR,NULL,ISTAT) IF (BLINK) THEN * *blink funny format: conversion to (+,-) microns XBLINK = XCUR YBLINK = YCUR CALL BLINKXY(XBLINK,YBLINK,XCUR,YCUR,1) 9999 Format(80a) ENDIF *--- Compute mean (x,y) for each object and R.M.S. in X,Y,R IF (ILOOP .LE. 1) THEN * *first (X,Y) AFOR = ACUR XTEMP(1) = XCUR YTEMP(1) = YCUR IX = 1 GO TO 62 ENDIF IF (ACUR .EQ. AFOR) THEN * *still same object IX = IX + 1 XTEMP(IX) = XCUR YTEMP(IX) = YCUR ENDIF IF (ACUR .NE. AFOR .OR. ILOOP .EQ. NROWD) THEN * *New object or last measurement c- *--> do means for the one before CALL MEAN(XTEMP,YTEMP,IX,XMEAN,YMEAN,NOUT,RESA,SIG) IF (IX.GT.2) THEN DO I = 1,3 RESB(I) = RESB(I) + RESA(I) ENDDO IRES = IRES + IX IRESNO = IRESNO + 1 ENDIF XTEMP(1) = XCUR YTEMP(1) = YCUR ITOT = ITOT + 1 * *Maximum # of objects is MAXOBJ IF (ITOT.GT.MAXOBJ) THEN WRITE(MSG,9074) CALL STTPUT(MSG,ISTAT) 9074 FORMAT('*error* too many objects ! ') CALL STSEPI STOP ENDIF PRXY(ITOT,1) = XMEAN/1.0D3 PRXY(ITOT,2) = YMEAN/1.0D3 * *conversion Text ident --> Nr ident IID = 0 628 IF (IID .LT. 16) THEN IID = IID+1 ELSE WRITE(MSG,*) '*warning* identiers must contain numbers' CALL STTPUT(MSG,ISTAT) WRITE(MSG,*) ' identifier set to 0' CALL STTPUT(MSG,ISTAT) NFOR = 0 GOTO 629 ENDIF READ(AFOR(IID:16),*,err=628) NFOR WRITE(MSG,9076) NFOR,XMEAN,YMEAN,IX+NOUT,RESA,NOUT CALL STTPUT(MSG,ISTAT) 9076 FORMAT(I8,2F12.3,I5,3F8.3,I5) 629 MNOST(ITOT) = NFOR AFOR = ACUR IX = 1 ENDIF * *go to read next line 62 CONTINUE * *file is finished 63 NOBJ = ITOT WRITE(MSG,9080) NOBJ CALL STTPUT(MSG,ISTAT) 9080 FORMAT(' No. of objects in (X,Y)-table: ',I5) * CALCULATE MEAN R.M.S. FOR MEASURED (X,Y) IN MICRONS IF (IRES.GT.0) THEN DO I = 1,3 RESB(I) = RESB(I)/FLOAT(IRESNO) ENDDO WRITE(MSG,9082) RESB CALL STTPUT(MSG,ISTAT) 9082 FORMAT(' Mean R.M.S. (microns or pix) in (X,Y,R) :',3f8.2) WRITE(MSG,9182) IRES,IRESNO CALL STTPUT(MSG,ISTAT) 9182 FORMAT(' based on ',i3,' (X,Y) measurements of ',i3, > ' objects ') ELSE WRITE(MSG,9083) CALL STTPUT(MSG,ISTAT) 9083 FORMAT('*info* Not enough measurements of same object to ', > ' calculate mean r.m.s. values ') ENDIF *---------------------------------------------------------------------- * * STANDARD STARS *---------------------------------------------------------------------- CALL STTPUT(' ',ISTAT) CALL STTPUT(' ',ISTAT) WRITE(MSG,9183) 9183 FORMAT(' --- STANDARD STARS ---') CALL STTPUT(MSG,ISTAT) NOSAO = 0 IF (.NOT. CENTER) THEN ALFA0 = 0. DELTA0 = 0. ENDIF DEP = EPOCH - EP1 C this loop from ILOOP = 1 to ILOOP = number of standards file DO 111 ILOOP = 1, NROWS * * Read the table: * *is the current star selected (= measured and non deleted) IF ( SELECT ) THEN CALL TBSGET(TIDS, ILOOP, ROWSEL, ISTAT) IF (.NOT. ROWSEL) GOTO 111 ENDIF * * Identifier (ppm #) and magnitude CALL TBERDI(TIDS,ILOOP,NCSPPM,NOS,NULL,ISTAT) CALL TBERDD(TIDS,ILOOP,NCSMAG,STM,NULL,ISTAT) * * R_A and Dec are BOTH read in degree CALL TBERDD(TIDS,ILOOP,NCSR_A,R_A,NULL,ISTAT) CALL TBERDD(TIDS,ILOOP,NCSDEC,DEC,NULL,ISTAT) * * pma and pmd are read in arcsec/year CALL TBERDD(TIDS,ILOOP,NCSPMA,PMA,NULL,ISTAT) CALL TBERDD(TIDS,ILOOP,NCSPMD,PMD,NULL,ISTAT) 9115 FORMAT(I6,F6.2,2I3,F7.3,A2,I2,I3,F6.2,2I6) * *cross check the measurements and the standards DO 120 I = 1,NOBJ C write(6,*) ILOOP,NOS,I,MNOST(i) IF (NOS.EQ.MNOST(I)) then * * the current object is a std! * * converts the coordinates 125 PRSTMAG(I) = STM PRST(I,1) = R_A*DOPI18 135 PRST(I,2) = DEC*DOPI18 * * and propoer motion * * PM_RA(radians) = pma(arcsec)/3600/cos(delta) * pi/180 PMAL = PMA*DOPI18 / DO3600/DCOS(PRST(I,2)) * * PM_Dec(radians) = pmd(arcsec) PMDEL = PMD*DOPI18/DO3600 * * coordinates a the epoch of the measurements PRST(I,1) = PRST(I,1)+PMAL*DEP PRST(I,2) = PRST(I,2)+PMDEL*DEP * ( Original pos1 was using plate center read from the standard stars * file header as reference for the alpha < or > 360 deg dilemma. * This version uses the 1st standard found ) IF (NOSAO.EQ.0) THEN * * First std: we keep it as ref. ALFA00 = PRST(I,1) ELSE * * adjust the alpha if no in the same region as the ref. IF(ALFA00.LT.PI/2D0.AND.PRST(I,1).GT.1.5D0*PI) > PRST(I,1) = PRST(I,1) - 2.0D0*PI IF(ALFA00.GT.1.5D0*PI.AND.PRST(I,1).LT.PI/2D0) > PRST(I,1) = PRST(I,1) + 2.D0*PI ENDIF IF (IFOUND(I).EQ.0)THEN * * this standard was not found before: * * so we have one more NOSAO = NOSAO + 1 IFOUND(I) = 1 IF(.NOT. CENTER)THEN ALFA0 = ALFA0 + PRST(I,1) DELTA0 = DELTA0 + PRST(I,2) ENDIF * * write in the table that the std was found CALL TBEWRI(TIDS,ILOOP,NCSSTD,1,ISTAT) * * store the row for further use ISTDR(I) = ILOOP ENDIF * * found, so it is not necessary to continue to search GO TO 111 ENDIF 120 CONTINUE 111 CONTINUE IF (NOSAO .EQ. 0) THEN CALL STTPUT('*ERROR* No standard star??',ISTAT) STOP ENDIF * *compute the plate/field center: IF (.NOT. CENTER) THEN * * Center is mean of the standards ALFA0 = ALFA0/DFLOAT(NOSAO) DELTA0 = DELTA0/DFLOAT(NOSAO) ENDIF * * ELSE: the coordinates have been read from the parameters. COSD0 = DCOS(DELTA0) SIND0 = DSIN(DELTA0) CALL ADCON(ALFA0,DELTA0,A0,D0) IF (CENTER) THEN WRITE(MSG,9091) 9091 FORMAT(' Plate/field center (read from the parameters):') ELSE WRITE(MSG,9591) 9591 FORMAT(' Plate/field center (MEAN of the standards coord.):') ENDIF CALL STTPUT(MSG,ISTAT) WRITE(MSG,9191) A0 CALL STTPUT(MSG,ISTAT) DSIG='+' IF(DELTA0.LT.0.) DSIG='-' WRITE(MSG,9291) DSIG,D0 CALL STTPUT(MSG,ISTAT) WRITE(MSG,9391) EPOCH CALL STTPUT(MSG,ISTAT) WRITE(MSG,9491) EP1 CALL STTPUT(MSG,ISTAT) 9191 FORMAT(' Alpha: ',F5.0,1X,F3.0,1X,F6.3) 9291 FORMAT(' Delta: ',1A,F4.0,1X,F3.0,1X,F5.2) 9391 FORMAT(' Plate epoch: ',F8.3) 9491 FORMAT(' Standard stars catalogue epoch: ',F8.3) * *decode the xy transformation terms selection flags 150 NX = 0 NY = 0 DO I = 1,9 READ(STRING_TERM(I:I),*)IDX(I) READ(STRING_TERM(I+10:I+10),*)IDY(I) IF (IDX(I).NE.0) NX=NX+1 IF (IDY(I).NE.0) NY=NY+1 ENDDO 165 FORMAT(1X,' CHECK ',9I2,1X,9I2) *--- Check if any stars deleted and delete it if yes 180 L = 0 ERXY(1) = DO0 ERXY(2) = DO0 ERXY(3) = DO0 DO 200 I = 1,NOBJ LOG = (IFOUND(I) .EQ. 1) IF (NKICK.NE.0) THEN DO 185 J = 1,NKICK IF (MNOST(I).EQ.KICK(J)) LOG = .FALSE. 185 CONTINUE ENDIF 190 IF ( LOG ) THEN * * this star is an undeleted standard L = L+ 1 NOST(L) = MNOST(I) XY(L,1) = PRXY(I,1) XY(L,2) = PRXY(I,2) ST(L,1) = PRST(I,1) ST(L,2) = PRST(I,2) STMAG(L) = PRSTMAG(I) ISTDR2(L) = ISTDR(I) 9195 FORMAT(2I7,6F15.10/) ENDIF 200 CONTINUE 205 NSTAN = NOSAO - NKICK Coli WRITE(MSG,9210) NOSAO,NKICK,NSTAN Coli CALL STTPUT(MSG,ISTAT) 9210 FORMAT(' Standard stars: ',I3,'. Manually removed: ', $ I2,'. Now left: ',I3) IF (NSTAN.NE.L) THEN WRITE(MSG,9215) NSTAN,L CALL STTPUT(MSG,ISTAT) ENDIF 9215 FORMAT('*ERROR* Wrong No. of standard stars !',2I6) IF (ITRY.GT.0) GO TO 325 C C OUTPUT STANDARD STAR FILE HEADER C CALL STTPUT(' ',ISTAT) CALL STTPUT(' ',ISTAT) WRITE(MSG,9315) CALL STTPUT(MSG,ISTAT) 9315 FORMAT(' --- ASTROMETRIC RESOLUTION ---') c$$$ WRITE(MSG,9090) NOFI c$$$ CALL STTPUT(MSG,ISTAT) c$$$ 9090 FORMAT(' FIELD NO.: ',I5) c$$$ WRITE(MSG,9190) IA3,A0(3) c$$$ CALL STTPUT(MSG,ISTAT) c$$$ 9190 FORMAT(' ALPHA: ',2I3,F7.3) c$$$ WRITE(MSG,9290) ID3,D0(3) c$$$ CALL STTPUT(MSG,ISTAT) c$$$ 9290 FORMAT(' DELTA: ',2I3,F7.2) c$$$ WRITE(MSG,9390) EPOCH c$$$ CALL STTPUT(MSG,ISTAT) c$$$ 9390 FORMAT(' PLATE EPOCH: ',F8.3) c$$$ WRITE(MSG,9490) EP1 c$$$ CALL STTPUT(MSG,ISTAT) c$$$ 9490 FORMAT(' STANDARD STARS CATALOGUE EPOCH: ',F8.3) 325 M = NSTAN DO 330 I = 1,NSTAN COSRH = COSD0*DCOS(ST(I,2))*DCOS(ST(I,1)-ALFA0)+ 1 SIND0*DSIN(ST(I,2)) SINRH = DSQRT(1.d0-COSRH**2) SINPH = DCOS(ST(I,2))*DSIN(ST(I,1)-ALFA0)/SINRH COSPH = (DSIN(ST(I,2))-SIND0*COSRH)/(COSD0*SINRH) * ( If plate curved (Schmidt) use distance along arc, otherwise use * linear distance in tangential plane ) IF (SCHMIDT) THEN RH = DATAN(SINRH/COSRH) ELSE RH = SINRH/COSRH ENDIF AKSI = RH*SINPH/DOPI18 ANU = RH*COSPH/DOPI18 XR = XY(I,1) YR = XY(I,2) AR(I,1) = XR AR(I,2) = YR AR(I,3) = AKSI AR(I,4) = ANU C write(6,*) i,(AR(I,io),io=1,4) 330 CONTINUE DO 650 IXY = 1,2 IF (IXY.EQ.2) then N = NY + 1 WRITE(MSG,9350) CALL STTPUT(' ',ISTAT) CALL STTPUT(MSG,ISTAT) 9350 FORMAT(' *- Y-Transformation') else N = NX + 1 WRITE(MSG,9340) CALL STTPUT(' ',ISTAT) CALL STTPUT(MSG,ISTAT) 9340 FORMAT(' *- X-Transformation') endif C C COMPUTE THE MATRICES C 355 L=N-1 SW=DBLE(FLOAT(M)) IZERO=32 IDSK=1 DO 360 I1=1,N SIGMA(I1)=DO0 DO 360 I1B=I1,N 360 SSX(I1,I1B)=DO0 DO 410 I3=1,M CALL RAR(I3,0,N,KK,IXY,IDX,IDY,X,AR) DO 380 I10=1,N SIGMA(I10)=SIGMA(I10)+X(I10) DO 380 I10B=I10,N 380 SSX(I10,I10B)=SSX(I10,I10B)+X(I10)*X(I10B) 410 CONTINUE WRITE(MSG,9420) M,N CALL STTPUT(MSG,ISTAT) 9420 FORMAT(' No of data = ',I3,' No of variables = ',I2) 430 RSW=DO1/SW DO 440 I4=1,N 440 XBAR(I4)=SIGMA(I4)*RSW DO 450 I5=1,N DO 450 I5B=I5,N 450 S(I5,I5B)=SSX(I5,I5B)-SW*XBAR(I5)*XBAR(I5B) 460 DO 470 I6=1,N SIY(I6)=S(I6,N) AH = S(I6,I6) 470 SIGMA(I6)=DSQRT(AH) DO 480 I11=1,L DO 480 I11B=I11,L SSINV(I11,I11B)=S(I11,I11B) 480 SSINV(I11B,I11)=SSINV(I11,I11B) CALL QKDIN(SSINV,L) DO 490 I12=1,L 490 SP(I12)=SSINV(I12,I12) 500 DO 510 I19=1,L B(I19)=DO0 DO 510 I19B=1,L 510 B(I19)=B(I19)+SIY(I19B)*SSINV(I19,I19B) CALL QKDIN(SSINV,L) 520 DO 530 I7=1,L K=I7+1 DO 530 I7B=K,N 530 R(I7,I7B)=S(I7,I7B)/(SIGMA(I7)*SIGMA(I7B)) K=N-2 540 DO 550 I13=1,L R(I13,I13)=DO1 DO 550 I13B=I13,L SSINV(I13,I13B)=R(I13,I13B) 550 SSINV(I13B,I13)=SSINV(I13,I13B) CALL QKDIN(SSINV,L) 560 SSR=DO0 SUM=DO0 DO 570 I14=1,L SUM=SUM+B(I14)*XBAR(I14) 570 SSR=SSR+B(I14)*SIY(I14) CONST =XBAR(N)-SUM SSDR=SIY(N) -SSR RSQD=SSR/SIY(N) ITDF=M-1 KDF=N-1 MK1DF=M-KDF-1 CHG=DBLE(FLOAT(KDF)) XMSDR=SSR/CHG CHG=DBLE(FLOAT(MK1DF)) XMSDD=SSDR/CHG ROFF=XMSDR/XMSDD VE=XMSDD IF (VE.LE.0.0) VE = 0.0 SDE=DSQRT(VE) DO 580 I15=1,L 580 VBI(I15)=VE*SP(I15) WRITE(MSG,9590) CONST CALL STTPUT(MSG,ISTAT) 9590 FORMAT(' Constant term = ',E15.8) WRITE(MSG,9600) CALL STTPUT(MSG,ISTAT) 9600 FORMAT(' Variable Coefficient' $ ,' Variation of coefficient') DO 595 J1 = 1,L I = KK(J1) WRITE(MSG,9601) (KTEXT(K,I),K=1,3), B(J1), VBI(J1) CALL STTPUT(MSG,ISTAT) 595 CONTINUE 9601 FORMAT(' Fct(',3A2,')',2X,E13.6,4x,E13.6) WRITE (MSG,9610) RSQD CALL STTPUT(MSG,ISTAT) 9610 FORMAT(' R square =',E18.5) IF (IXY.EQ.1) THEN SCALE = B(1)*3.6D3 ELSE SCALE = B(2)*3.6D3 ENDIF SCALE = DABS(SCALE) WRITE(MSG,9615) SCALE CALL STTPUT(MSG,ISTAT) 9615 FORMAT(' Approx. plate scale: ', F10.3, $ ' arcsec per 1000 Pixel or Micron ') IDSK=1 DO 630 I16=1,M V=DO0 CALL RAR(I16,0,N,KK,IXY,IDX,IDY,X,AR) DO 620 I17=1,L 620 V=V+B(I17)*X(I17) YHAT=CONST+V XY(I16,IXY) = YHAT DEV=X(N)-YHAT 630 CONTINUE IF (IXY.EQ.1) CX = CONST IF (IXY.EQ.2) CY = CONST DO 640 J=1,9 IF (IXY.EQ.1) BX(J) = B(J) IF (IXY.EQ.2) BY(J) = B(J) 640 CONTINUE 650 CONTINUE c WRITE(MSG,9660)NOFI,EPOCH c CALL STTPUT(MSG,ISTAT) 9660 FORMAT('***pipo*** FIELD NO.',I5,' EPOCH',F10.3, > ' Epoch of std stars cat') c WRITE(MSG,9760)A0,D0 c CALL STTPUT(MSG,ISTAT) 9760 FORMAT('***pipo*** FIELD CENTER ',3F4.0,1X,3F4.0) c WRITE(MSG,9860)NSTAN c CALL STTPUT(MSG,ISTAT) 9860 FORMAT('***pipo*** NO OF STANDARD STARS MEASURED',I5) CALL STTPUT(' ',ISTAT) CALL STTPUT(' ',ISTAT) WRITE(MSG,9960) CALL STTPUT(MSG,ISTAT) 9960 FORMAT('Row No Mag Measured' $ ' Catalogue Deviation' ) WRITE(MSG,9963) CALL STTPUT(MSG,ISTAT) 9963 FORMAT('<-> <----> <--> <---RA---> <---Dec---> ' $ ,'<---RA---> <---Dec---> -X-arcsec-Y-') cc WRITE(MSG,9661) cc CALL STTPUT(MSG,ISTAT) cc 9661 FORMAT(2(' Row No Mag ',4X,'(X)RMS(Y) ')) cc WRITE(MSG,9761) cc CALL STTPUT(MSG,ISTAT) cc 9761 FORMAT(2(23X,'Arcsec ')) SUMD1 = DO0 SUMD2 = DO0 DO 700 I = 1,NSTAN AKSI = XY(I,1) ANU = XY(I,2) CALL XYAD(AKSI,ANU,IA1,ID1,DA1S,DD1S,ITEGN1,ALFAS,DELTS, * COSD0,SIND0,ALFA0) CALL ADCON(ST(I,1),ST(I,2),ALFA,DELTA) DO 680 J=1,2 IA2(J) = IDINT(ALFA(J)+1.0D-3) 680 ID2(J) = IDINT(DELTA(J)+DSIGN(1.0D-3,DELTA(J))) ID2(1) = IABS(ID2(1)) ITEGN2 = ' +' IF (ST(I,2).LT.DO0) ITEGN2 = ' -' SUMD3 = (ALFAS - ST(I,1))*DCOS(ST(I,2))*DO3600/DOPI18 SUMD4 = (DELTS-ST(I,2))*DO3600/DOPI18 SUMD5 = DSQRT(SUMD3**2+SUMD4**2) SUMD1 = SUMD1 + SUMD3**2 SUMD2 = SUMD2 + SUMD4**2 * *write the errors in the standard table: CALL TBEWRD(TIDS,ISTDR2(I),NCSXER,SUMD3,ISTAT) CALL TBEWRD(TIDS,ISTDR2(I),NCSYER,SUMD4,ISTAT) IF (SUMD5.GT.ERXY(3)) THEN IERROW = ISTDR2(I) IERXY = NOST(I) ERXY(1) = SUMD3 ERXY(2) = SUMD4 ERXY(3) = SUMD5 endif 685 CONTINUE WRITE(MSG,9690)ISTDR2(I),NOST(I),STMAG(I), $ IA1(1),IA1(2),DA1S,ITEGN1, > ID1(1),ID1(2),DD1S,IA2(1),IA2(2),ALFA(3),ITEGN2,ID2(1), > ID2(2),DELTA(3),SUMD3,SUMD4 9690 FORMAT(I3,I7,F5.1,2(I3,I2,F6.3,A2,I2,I2,F6.2),2F7.3) CALL STTPUT(MSG,ISTAT) c$$$ WRITE(MSG,9691) ISTDR2(I),NOST(I),STMAG(I),SUMD3,SUMD4 c$$$ CALL STTPUT(MSG,ISTAT) c$$$ 9691 FORMAT(I5,I7,F4.1,2X,2F7.3,' ',$) c$$$C 9691 FORMAT(I6,F4.1,2X,2F7.3,' ',$) c$$$ IF (MOD(I,2).EQ.0) WRITE(MSG,9692) c$$$ CALL STTPUT(MSG,ISTAT) c$$$ 9692 FORMAT($) 700 CONTINUE RH = DBLE(FLOAT(NSTAN-1)) SUMD1 = DSQRT(SUMD1/RH) SUMD2 = DSQRT(SUMD2/RH) c$$$ WRITE(MSG,9692) c$$$ CALL STTPUT(MSG,ISTAT) WRITE(MSG,9710)SUMD1 CALL STTPUT(MSG,ISTAT) WRITE(MSG,9711)SUMD2 CALL STTPUT(MSG,ISTAT) 9710 FORMAT(' RMS in X-direction (arcsec):',F7.3) 9711 FORMAT(' RMS in Y-direction (arcsec):',F7.3) 740 IF (DABS(ERXY(1)).ge.DO1 .or. DABS(ERXY(2)).ge.DO1) THEN WRITE(MSG,742) IERXY,IERROW,ERXY(1),ERXY(2) CALL STTPUT(MSG,ISTAT) ENDIF 742 FORMAT(' *TIP* star No. ',i8,' (row ',i5,') has large RMS: ' > ,2f10.2) *--- write the transformation matrix in output keyword CALL STKWRD('BX',BX,1,9,IDUM,ISTAT) CALL STKWRD('BY',BY,1,9,IDUM,ISTAT) CALL STKWRD('CX',CX,1,1,IDUM,ISTAT) CALL STKWRD('CY',CY,1,1,IDUM,ISTAT) CALL STKWRD('AL_DE0',ALFA0,1,1,IDUM,ISTAT) CALL STKWRD('AL_DE0',DELTA0,2,1,IDUM,ISTAT) CALL STKWRD('AL_DE0',COSD0,3,1,IDUM,ISTAT) CALL STKWRD('AL_DE0',SIND0,4,1,IDUM,ISTAT) CALL STKWRI('NX',NX,1,1,IDUM,ISTAT) CALL STKWRI('NY',NY,1,1,IDUM,ISTAT) CALL STSEPI 9045 FORMAT(A50) 9050 FORMAT(2X,A50) END