C @(#)astrm.for 17.1.1.1 (ES0-DMD) 01/25/02 17:18:59 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 Massachusetss Ave, Cambridge, C MA 02139, USA. C C Corresponding 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 C +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ C.COPYRIGHT: Copyright (c) 1987 European Southern Observatory, C all rights reserved C.IDENTIFICATION: ASTRM C.PURPOSE: Routine for main program star written by P.B.Stetson C performs the gaussian astrometric fitting C.LANGUAGE: F77+ESOext C.AUTHOR: P. Stetsan C.ALGORITHM: see Stetson, P.B., 1979 Astron. J. 84 1149 C.INPUT/OUTPUT: -> KRX input array C -> JX low side local minimun C -> LX high side local minimum C <-> XI input/ouput estimates of the image centre C <- S1 error estimate of the output image centre C 0. fitting procedure failed or if C KOUNT .NE. 0 then KOUNT too small C <-> SX input/output estimate of the image width C <- DO output estimate of the image height C <- B output estimate of the background marginal iensity. C <-> KOUNT number of iterations required/ status C 0 fitting procedure failed C.COMMENTS: If the solution fails to converge, XI will come back with C the input value.... with S1 = 0. just in case you want to try C doing the photometry with the estimated centre from SERCH. C.KEYWORDS : photometry,astrometry C.VERSION: 8212?? P. Stetsan C.VERSION: 871123 R.H. Warmels ESO Fortran Conversion C.VERSION: ?????? H. Waldhausen Vax implementation C --------------------------------------------------------------------- SUBROUTINE ASTRM(KRX,JX,LX,XI,S1,SX,DO,B,KOUNT) C IMPLICIT NONE DOUBLE PRECISION KRX(1) INTEGER JX INTEGER LX REAL XI REAL S1 REAL SX REAL DO REAL B INTEGER KOUNT C INTEGER I INTEGER J INTEGER KNT INTEGER L REAL D, DX, DELX REAL EX REAL F, FF, FOLD REAL S, S2, SXSI REAL XO REAL Y DOUBLE PRECISION CTF(5) DOUBLE PRECISION CTC(5,5) DOUBLE PRECISION V(5) DOUBLE PRECISION DABS REAL C(3) C C C.. this part of the program uses the estimates of the image centre C.. and width from SERCH to derive initial estimates of the image C.. height and background density by linear least-squares. C IF ((LX-JX).GT.4) THEN XO = XI DELX = LX - JX + 1 XO = XI DELX = LX - JX + 1 S1 = 0.0 S2 = 0.0 DO = 0.0 B = 0.0 SXSI = 0.5/(SX*SX) C DO 10 I = JX,LX DX = I - XO EX = EXP(-DX*DX*SXSI) DO = DO + EX*KRX(I) B = B + EX*EX S1 = S1 + KRX(I) S2 = S2 + EX 10 CONTINUE C F = LX - JX + 1 DO = (DO-S1*S2/F)/ (B-S2*S2/F) B = (S1-DO*S2)/F IF (DO.LT.20.0) D = 20.0 IF (SX.LT.2.0) SX = 2.0 FOLD = 1.0E38 C C.. main iteration loop DO 90 KNT = 1,KOUNT ! iteration loop S1 = 0.0 S2 = 0.0 FF = 0.0 SXSI = 1.0/ (SX*SX) DO 30 I = 1,4 CTF(I) = 0.0D0 DO 20 J = I,4 CTC(I,J) = 0.0D0 20 CONTINUE 30 CONTINUE C C.. set normal equations DO 60 L = JX,LX ! equation loop DX = L - XO EX = 0.5*DX*DX*SXSI C(2) = EXP(-EX) S = DO*C(2) Y = S + B F = Y - KRX(L) C(1) = DX*SXSI*S EX = C(1)*C(1) S1 = S1 + EX*F*F S2 = S2 + EX FF = FF + F*F C(3) = C(1)*DX/SX DO 50 I = 1,3 CTF(I) = CTF(I) + F*C(I) CTC(I,4) = CTC(I,4) + C(I) DO 40 J = I,3 CTC(I,J) = CTC(I,J) + C(I)*C(J) 40 CONTINUE 50 CONTINUE CTF(4) = CTF(4) + F ! equation loop 60 CONTINUE C CTC(4,4) = FLOAT(LX-JX+1) DO 80 I = 2,4 L = I - 1 DO 70 J = 1,L CTC(I,J) = CTC(J,I) 70 CONTINUE 80 CONTINUE C C.. solve normal equations CALL INVRS(CTC,4) CALL VMUL(CTC,4,CTF,V) C C.. if DA is the correction to the estimate of A, then replace DA by C.. DA/(1+(DA/X)) = DA' prevents an ill-determined solution from C.. jumping around too much in subsequent iterations. DA' is less C.. than X, and DA' goes to DA for DA much less than X. C.. If X = A/R for some R greater than one, then an intrinsically C.. positive parameter will never be permitted to pass through zero. C.. This can be helpful for parameters like image width or height. C XO = XO - V(1)/ (1.0+DABS(V(1)/5.0)) DO = DO - V(2)/ (1.0+2.0*DABS(V(2)/DO)) SX = SX - V(3)/ (1.0+2.0*DABS(V(3)/SX)) B = B - V(4)/ (1.0+2.0*DABS(V(4)/B)) C C.. failure IF ((DO.LT.20.0) .OR. (SX.LT.1.0) .OR. + (SX.GT.0.5*DELX)) THEN KOUNT = 0 S1 = 0.0 RETURN END IF C C.. success IF ((DABS(V(1)).LT.0.01) .AND. + ((FOLD-FF)/FOLD.LT.0.05)) THEN XI = XO S1 = S1*CTC(1,1)*DELX/ (S2* (DELX-4.0)) RETURN END IF C FOLD = FF 90 CONTINUE ! iteration loop C C.. number of iterations too small S1 = 0.0 C C.. undersampled ELSE KOUNT = 0 S1 = 0.0 END IF C.. exit RETURN END