C @(#)sgfpi.for 17.1.1.1 (ES0-DMD) 01/25/02 17:19:21 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 C.IDENTIFICATION: C PROGRAM SGFPI VERSION 1.4 FEB 9, 1988 C P.GROSBOL ESO - GARCHING C C.PURPOSE: C PROGRAM FITS THE POSITION AND INCLINATION ANGLE OF A GALAXY SO C THAT THE PHASES OF 2ND AND 4TH HARMONIC HAVE THE MINIMUM SIGMA. C C.INPUT/OUTPUT: C THE PROGRAM RUNS UNDER MIDAS AND USES THE FOLLOWING KEYWORDS: C IN_A : NAME OF IMAGE FILE (I) C POSI(4) : REAL INPUT PARAMETERS DEFINING GALAXY POSITION (I/O) C 1:X CENTER POSITION, 2:Y CENTER POSITION, C 3:POSITON ANGLE (DEG), 4:INCLINATION ANGLE (DEG) C RANGE(3) : REAL INPUT PARAMETERS DEFINING RADIAL RANGE (I) C 1:INNER RADIUS, 2:OUTER RADIUS, 3:RADIAL STEP C REGION(2) : THE REGION OF PA AND I TO SEARCH FOR A MINIMUM (I) C 1:PA REGION (DEG), 2:I REGION (DEG) C C--------------------------------------------------------------------- PROGRAM SGFPI C IMPLICIT NONE INTEGER MR INTEGER MS INTEGER MFC PARAMETER (MR=100) PARAMETER (MS=7) PARAMETER (MFC=5) C INTEGER NPIX(3),PNTR,IMNO,KUNIT(4),INULL INTEGER MADRID(1) INTEGER IERR, IAV, IP, II, I INTEGER IX, IDF INTEGER MC, NR, NA C REAL AMP(MFC),PHA(MFC),PH(MFC), , RANGE(3),POSI(4),REG(2), , RL(MR),P2(MR),P4(MR),W(MR), , A(MS,MS),V(MS,MS),S(MS,MS),DA(MS,MS),DB(MS,MS), , SA(MS,MS),SB(MS,MS),RR(MS,MS),DIA(MS) REAL PIQ, PI, PAS, PA, PAR, PI2, PIH REAL DTR, DR, DPA, DAI, DIF REAL RI, RO, RM, R, RP, RD REAL XC, YC REAL AIS, AIR REAL A2, A4 REAL AN, AI, AD REAL B2, B4 REAL R2, R4 REAL BD REAL SM REAL SA2, SB2, SA4, SB4 REAL SAD, SBD, SS2, SS4 REAL XM, WM C DOUBLE PRECISION START(3),STEP(3) C LOGICAL LFPHA C CHARACTER*60 INAME CHARACTER*64 UNIT CHARACTER*72 IDEN CHARACTER*80 MSG C INCLUDE 'MID_INCLUDE:ST_DEF.INC' COMMON /VMR/MADRID INCLUDE 'MID_INCLUDE:ST_DAT.INC' C C INITIATE MIDAS C CALL STSPRO('SGFPI') C C READ PARAMETERS C CALL STKRDC('IN_A',1,1,60,IAV,INAME,KUNIT,INULL,IERR) CALL STKRDR('INPUTR',5,3,IAV,RANGE,KUNIT,INULL,IERR) CALL STKRDR('INPUTR',1,4,IAV,POSI,KUNIT,INULL,IERR) CALL STKRDR('INPUTR',8,2,IAV,REG,KUNIT,INULL,IERR) C C OPEN IMAGE FRAME C CALL STIGET(INAME,D_R4_FORMAT,F_I_MODE,1,3,NA,NPIX,START, 2 STEP,IDEN,UNIT,PNTR,IMNO,IERR) IF (NA.NE.2 .OR. IERR.NE.0 .OR. . ABS(SNGL(STEP(1))/SNGL(STEP(2))-1.0).GT.1E-5) THEN CALL STTPUT('*** FATAL: Wrong image format',IERR) CALL STSEPI ENDIF C C CHECK PARAMETERS AND CONVERT TO PIXEL VALUES C LFPHA = .TRUE. PIQ = ATAN(1.0) PIH = 2.0 * PIQ PI = 4.0 * PIQ PI2 = 8.0 * PIQ DTR = PI / 180.0 RI = MAX( 0.0,ABS(RANGE(1)) ) RO = MAX( RI,ABS(RANGE(2)) ) RM = MIN( ABS(POSI(1)-SNGL(START(1))), , ABS(POSI(2)-SNGL(START(2))) ) RM = MIN( RM,ABS(POSI(1)-SNGL(START(1))-NPIX(1)*SNGL(STEP(1))) ) RM = MIN( RM,ABS(POSI(2)-SNGL(START(2))-NPIX(2)*SNGL(STEP(2))) ) RO = MIN( RM-2.0*ABS(SNGL(STEP(1))),RO ) DR = ABS( RANGE(3) ) XC = (POSI(1)-SNGL(START(1)))/SNGL(STEP(1)) + 1.0 YC = (POSI(2)-SNGL(START(2)))/SNGL(STEP(2)) + 1.0 NR = MIN( MR,INT((RO-RI)/DR) ) DPA = REG(1) / (MS-1) DAI = REG(2) / (MS-1) PAS = POSI(3) - 0.5*REG(1) AIS = POSI(4) - 0.5*REG(2) DO 60, I = 1,MS DIA(I) = AIS + (I-1)*DAI 60 CONTINUE C C GO THROUGH ALL POINTS IN THE MESH C PA = PAS DO 1000, IP = 1,MS AI = AIS DO 1100, II = 1,MS C C GET FOURIER COEFFICIENTS FROM EACH RADIUS C XM = 0.0 SM = 0.0 WM = 0.0 R = RI PAR = DTR * PA AIR = DTR * AI DO 1200, I = 1,NR RP = R / SNGL(STEP(1)) RL(I) = LOG(R) P2(I) = 0.0 P4(I) = 0.0 W(I) = 0.0 CALL APFFTC(MADRID(PNTR),NPIX(1),NPIX(2),XC,YC, , RP,PAR,AIR,MFC,MC,AMP,PHA,AN) IF (MC.GE.MFC) THEN IF (LFPHA) THEN LFPHA = .FALSE. PH(1) = 0.0 DO 1250, IX = 2, MFC DIF = PI2 / FLOAT(IX-1) IDF = INT((PHA(IX)-PH(IX-1))/DIF+100.5)-100 PH(IX) = PHA(IX) - FLOAT(IDF)*DIF 1250 CONTINUE ENDIF CALL PHAMOD(PHA,PH,MFC) P2(I) = PHA(3) P4(I) = PHA(5) W(I) = 1.0 DIF = ATAN(TAN(2.0*(P2(I)-P4(I)))) / DTR XM = XM + DIF SM = SM + DIF*DIF WM = WM + 1.0 ENDIF R = R + DR 1200 CONTINUE C C COMPUTE MEAN AND RMS C CALL LFIT(RL,P2,W,NR,1,A2,SA2,B2,SB2,R2) CALL LFIT(RL,P4,W,NR,1,A4,SA4,B4,SB4,R4) CALL LFIT(P2,P4,W,NR,1,AD,SAD,BD,SBD,RD) XM = XM / WM SM = SQRT( MAX(0.0,SM/WM - XM*XM) ) A(IP,II) = XM S(IP,II) = SM DA(IP,II) = AD DB(IP,II) = BD V(IP,II) = RD SS2 = 0.0 SS4 = 0.0 DO 1300, I = 1,NR SS2 = SS2 + (P2(I)-A2-B2*RL(I))**2 SS4 = SS4 + (P4(I)-A4-B4*RL(I))**2 1300 CONTINUE SA(IP,II) = SQRT(SS2/NR) SB(IP,II) = SQRT(SS4/NR) RR(IP,II) = R2 AI = AI + DAI 1100 CONTINUE PA = PA + DPA 1000 CONTINUE C C ESTIMATE POSITION OF MINIMUM C WRITE(MSG,620) DIA 620 FORMAT(' INCL. : ',5X,7F8.2) CALL STTPUT(MSG,IERR) PA = PAS CALL STTPUT('REGR. COEF. P2-P4',IERR) DO 2000, IP = 1,MS WRITE(MSG,640) PA,(V(IP,II),II=1,MS) 640 FORMAT(' PA : ',F7.1,3X,7F8.4) 642 FORMAT(' PA : ',F7.1,3X,7F8.2) CALL STTPUT(MSG,IERR) PA = PA + DPA 2000 CONTINUE CALL STTPUT('SIGMA(DIF.)',IERR) PA = PAS DO 2100, IP = 1,MS WRITE(MSG,642) PA,(S(IP,II),II=1,MS) CALL STTPUT(MSG,IERR) PA = PA + DPA 2100 CONTINUE CALL STTPUT('MEAN(DIF.)',IERR) PA = PAS DO 2200, IP = 1,MS WRITE(MSG,642) PA,(A(IP,II),II=1,MS) CALL STTPUT(MSG,IERR) PA = PA + DPA 2200 CONTINUE CALL STTPUT('B(P2-P4)',IERR) PA = PAS DO 2300, IP = 1,MS WRITE(MSG,640) PA,(DB(IP,II),II=1,MS) CALL STTPUT(MSG,IERR) PA = PA + DPA 2300 CONTINUE CALL STTPUT('SIGMA P2',IERR) PA = PAS DO 2400, IP = 1,MS WRITE(MSG,640) PA,(SA(IP,II),II=1,MS) CALL STTPUT(MSG,IERR) PA = PA + DPA 2400 CONTINUE CALL STTPUT('SIGMA P4',IERR) PA = PAS DO 2500, IP = 1,MS WRITE(MSG,640) PA,(SB(IP,II),II=1,MS) CALL STTPUT(MSG,IERR) PA = PA + DPA 2500 CONTINUE CALL STTPUT('REGR. COEF. R2',IERR) PA = PAS DO 2600, IP = 1,MS WRITE(MSG,640) PA,(RR(IP,II),II=1,MS) CALL STTPUT(MSG,IERR) PA = PA + DPA 2600 CONTINUE C C FINISHED - EXIT C CALL STSEPI END