C psf 96,12 PARAMETER (MAXEXP=10, MAXPAR=6) PARAMETER (NOPT=20, MAXPSF=207) real opt(nopt) real psf(MAXPSF,MAXPSF,MAXEXP) real par(MAXPAR) integer ifail,iflag(MAXPSF) real a_b(4096*4032) character*60 f1,f2,f3,ch50 character*80 head(72) logical logi if(iargc().lt.1)then write(*,*) write(*,*) write(*,*)' ********** DAO_PSF 2010.3 ***********' write(*,*) write(*,*) write(*,*)' Usage: u_psf4k file.fit [-1,0,1,2] [!]' write(*,*) write(*,*) write(*,*)' if add !, u_psf4k DONOT delete any star' stop endif call getarg(1,f1) 10 k=index(f1,'.') if(k.eq.0)f1=f1(1:lnblnk(f1))//'.fit' write(*,'(a)')f1 inquire(file=f1,exist=logi) if(.not.logi)stop ' file not found !' f3='bok.par' inquire(file=f3,exist=logi) if(.not.logi)f3='/vega2/rhbin/bok.par' open(1,file=f3,status='old') do 20 i=1,nopt 21 read(1,'(a)')f3 if(f3(1:2).eq."**")goto 21 20 read(f3(1:),'(29x,f9.2)')opt(i) c write(*,'(7f10.2)')opt close(1) open(1,file=f1,status='old',access='direct',recl=2880*2) read(1,rec=1)head ipos=indexpos(head,'SEEING ') ! SEEING if(ipos.lt.72)then read(head(ipos)(25:),*)x ! 2001,10 if(x.ne.0)opt(5)=x/0.45*2 write(*,*)'FWHM(seeing in pixel)=',opt(5) endif close(1) call readfits(f1,a_b,ncol,nrow) k=index(f1,'.') f2=f1(1:k)//'lst' open(1,file=f2,status='old') ! 1 -- *.lst f2=f1(1:k)//'ap' open(2,file=f2,status='old') ! 2 -- *.ap f2=f1(1:k)//'psf' open(3,file=f2,status='UNKNOWN') ! 3 -- *.psf f2=f1(1:k)//'nei' open(4,file=f2,status='UNKNOWN') ! 4 -- *.nei if(iargc().gt.1)then call getarg(2,f2) read(f2(1:),*)opt(14) endif write(*,*)'variable PSF: ',opt(14) ifail=0 if(iargc().gt.2)then call getarg(3,f2) if(f2(1:1).eq.'!')ifail=-1 endif do i=1,MAXPSF iflag(i)=0 enddo call GETPSF (a_b,NCOL,NROW,PAR,PSF,OPT,NOPT,ifail,iflag) if(ifail.eq.1)then write(*,*)' *********** delete Bad PSF_star, redo !' k=index(f1,'.') f2=f1(1:k)//'lst' open(1,file=f2,status='old') ! 1 -- *.lst f3=f1(1:k)//'tmp' open(2,file=f3,status='UNKNOWN') i=0 40 i=i+1 read(1,'(a)',end=50)ch50 if(iflag(i).eq.0) write(2,'(a)')ch50 goto 40 50 close(1) close(2) call unlink(f2) call system("mv "//f3//" "//f2) goto 10 endif end SUBROUTINE GETPSF (PIC,NCOL,NROW,PAR,PSF,OPT,NOPT,ifail,iflag) IMPLICIT NONE C INTEGER MAXBOX, MAXPSF, MAXSTR, MAXN, MAXPAR, MAXEXP, NCOL, NROW INTEGER NOPT PARAMETER (MAXPAR=6, MAXPSF=207, MAXEXP=10) PARAMETER (MAXBOX=69, MAXSTR=50000, MAXN=200) integer ifail,iflag(MAXPSF) C CHARACTER LINE*80, LABEL*8 CHARACTER FLAG(MAXN)*1, ANSWER*1 DOUBLE PRECISION CC(2,2), VV(2), AA(2) REAL C(MAXEXP,MAXEXP), V(MAXEXP), A(MAXEXP) REAL CON(MAXPSF,MAXPSF), TERM(MAXEXP) REAL PIC(NCOL,NROW) REAL PAR(MAXPAR) REAL PSF(MAXPSF,MAXPSF,MAXEXP) REAL XCEN(MAXSTR), YCEN(MAXSTR), APMAG(MAXSTR), SKY(MAXSTR) REAL HPSF(MAXN), WEIGHT(MAXN), OPT(NOPT) REAL HJNK(MAXN), XJNK(MAXN), YJNK(MAXN) REAL PROFIL, BICUBC, SQRT INTEGER ID(MAXSTR), NTAB(-1:3), NPARAM LOGICAL IN(MAXPSF,MAXPSF), EDGE(MAXPSF,MAXPSF), SATR8D(MAXN) C DOUBLE PRECISION PSFMAG, DLOG10, DBLE REAL LOBAD, HIBAD, THRESH, AP1, PHPADU, READNS, RONOIS, . FWHM, WATCH, FITRAD, RADSQ, PSFRAD, PSFRSQ, FMAX, . DFDX, DFDY, SIG, DX, DY, RDX, RDY, XMID, YMID, DP, OLD, W, . WT, RADIUS, SCALE, SUMSQ, SUMN, DATUM, RSQ, DX2, DY2 REAL ERRMAG, CHI, SHARP, PERR, PKERR INTEGER I, J, K, L, N, LX, LY, MX, MY, NTOT, NSTAR, ISTAR, NL INTEGER IEXPAND, IFRAC, IPSTYP, NPSF, NPAR, ISTAT, IRADSQ, . NEXP, MIDDLE, MIDSQ, IDX, JDY, JDYSQ, IX, JY, ITER, NPASS INTEGER II, JJ, IR, JRADSQ, NPSFSQ LOGICAL STAROK, SATUR8, REDO logical logi C COMMON /ERROR/ PHPADU, RONOIS, PERR, PKERR C DATA NTAB / 0, 1, 3, 6, 10 / C C----------------------------------------------------------------------- C C SECTION 1 C C Set up the necessary variables, open the necessary files, read in the C relevant data for all stars. C FWHM = OPT(5) WATCH = OPT(11) FITRAD = OPT(12) RADSQ = FITRAD**2 PSFRAD = OPT(13) PSFRSQ = (PSFRAD+2.)**2 IEXPAND = NINT(OPT(14)) IFRAC = NINT(OPT(15)) NPASS = NINT(OPT(17)) C NPSF = 2*(NINT(2.*PSFRAD)+1)+1 PSFRAD = (REAL(NPSF-1)/2. - 1.)/2. IRADSQ = NINT(2.*PSFRAD)**2 JRADSQ = NINT(1.4*PSFRAD)**2 C C Ascertain the name of the input (aperture, usually) photometry file, C and read in the relevant data for all stars. C read(2,'(a)')i ! *.ap read(2,*) NL,I,I,LOBAD,HIBAD,THRESH,AP1,PHPADU,READNS,DX RONOIS = READNS**2 C I = 0 1010 I = I+1 IF (I .GT. MAXSTR) GO TO 1100 1020 CALL RDSTAR (2, NL, ID(I), XCEN(I), YCEN(I), APMAG(I), SKY(I)) IF (ID(I) .LT. 0) GO TO 1100 ! End-of-file encountered IF (ID(I) .EQ. 0) GO TO 1020 ! Blank line encountered GO TO 1010 C 1100 NTOT = I-1 close(2) C C----------------------------------------------------------------------- C C SECTIONS 2 C C Learn the name of each PSF star and find it in the star list. Display C a subarray around it if desired, and check it for invalid pixels. C NSTAR = 0 N = NPARAM(1, FWHM, LABEL, PAR, MAXPAR) C C NSTAR will point at the last PSF star in the stack after the PSF C stars have been brought to the top. C c read(1,'(a)')i c. read(1,'(a)')i c. read(1,'(a)')i 2000 CONTINUE read(1,*,end=3000,err=2000) I, DX, DY, RDX, RDY ! read *.lst IF (I .EQ. 0) GO TO 2000 C DO 2010 ISTAR=1,NTOT IF (ID(ISTAR) .EQ. I) GO TO 2020 2010 CONTINUE WRITE (6,701) I, ' was not found.' 701 FORMAT (1X, I5, A) GO TO 2000 C C If a given star appears in the .LST file more than once, employ it C only once. C 2020 IF (ISTAR .LE. NSTAR) GO TO 2000 IF ( ( INT(XCEN(ISTAR)-FITRAD) .LT. 0 ) .OR. . ( INT(XCEN(ISTAR)+FITRAD) .GT. NCOL ) .OR. . ( INT(YCEN(ISTAR)-FITRAD) .LT. 0 ) .OR. . ( INT(YCEN(ISTAR)+FITRAD) .GT. NROW ) ) THEN WRITE (6,701) ID(ISTAR), ' is too near edge of frame.' GO TO 2000 END IF C C Define the subarray containing the star and search it for "bad" C pixels. This will include a three-pixel border around the PSF box C itself, to guarantee valid interpolations to within the PSF box, and C to allow for a possible one-pixel shift in the star's centroid. C LX = MAX0( 1, INT(XCEN(ISTAR)-PSFRAD)-2 ) LY = MAX0( 1, INT(YCEN(ISTAR)-PSFRAD)-2 ) MX = MIN0(NCOL, INT(XCEN(ISTAR)+PSFRAD)+3 ) MY = MIN0(NROW, INT(YCEN(ISTAR)+PSFRAD)+3 ) C FMAX=-32768. STAROK=.TRUE. SATUR8=.FALSE. C 2040 REDO = .FALSE. DO 2070 J=LY,MY DP = (REAL(J)-YCEN(ISTAR))**2 DO 2050 I=LX,MX DX = REAL(I)-XCEN(ISTAR) RSQ = DX**2 + DP IF (RSQ .GT. PSFRSQ) THEN IF (DX .GE. 0.) GO TO 2070 GO TO 2050 END IF DATUM = PIC(I,J) IF ((DATUM .LT. LOBAD) .OR. (DATUM .GT. HIBAD)) THEN C C There is a defective pixel. C IF (RSQ .LE. RADSQ) THEN C C It is inside the fitting radius. C IF (DATUM .LT. LOBAD) THEN C C The pixel is bad low, so reject the star. C WRITE (6,701) ID(ISTAR), . ' is not a good star.' write(*,*)datum,lobad GO TO 2000 ELSE IF (.NOT. SATUR8) THEN C C The pixel is bad high, so presume the star is saturated. C WRITE (6,701) ID(ISTAR), . ' is saturated.' IF (WATCH .GT. -0.5) THEN write(*,'(a,$)')'Use it anyway? ' read(*,'(a)')ANSWER IF (ANSWER .EQ. 'Y') THEN SATUR8 = .TRUE. ELSE GO TO 2000 END IF ELSE SATUR8 = .TRUE. END IF END IF GO TO 2050 END IF C C The bad pixel isn't inside the fitting radius. C IF (STAROK) THEN STAROK = .FALSE. IF (WATCH .GT. -0.5) THEN WRITE (6,628) ID(ISTAR), DATUM, 7, I, J 628 FORMAT (1X, I5, ' has a bad pixel:', F9.1, . ' at position', A1, 2I5) END IF ELSE IF (WATCH .GT. -0.5) . WRITE (6,628) ID(ISTAR), DATUM, 0, I, J END IF C C Replace the defective pixel with an average of surrounding C pixels. C DATUM = 0. SUMN = 0. DO L=MAX0(1,J-1),MIN0(NROW,J+1) DO K=MAX0(1,I-1),MIN0(NCOL,I+1) IF ((PIC(K,L) .GE. LOBAD) .AND. . (PIC(K,L) .LE. HIBAD)) THEN DATUM = DATUM + PIC(K,L) SUMN = SUMN+1. END IF END DO END DO C C If there are fewer than three valid neighboring pixels, do not C fudge the pixel yet. Instead, we will go through again (and again C if necessary) and allow the fudge to eat into the defect from the C outside. C IF (SUMN .LT. 2.5) THEN REDO = .TRUE. ELSE PIC(I,J) = DATUM/SUMN END IF ELSE C C The pixel wasn't bad. C IF (DATUM .GT. FMAX) FMAX=DATUM END IF 2050 CONTINUE 2070 CONTINUE IF (REDO) GO TO 2040 C IF (STAROK) THEN IF ((.NOT. SATUR8) .AND. (WATCH .LT. -0.5) .AND. . (WATCH .GT. -1.5)) THEN C WRITE (6,701) ID(ISTAR), ' seems fine.' END IF ELSE IF (.NOT. SATUR8) THEN IF ((WATCH .LT. -0.5) .AND. (WATCH .GT. -1.5)) THEN WRITE (6,701) ID(ISTAR), ' has bad pixels.' ELSE IF (WATCH .GT. -0.5) THEN write(*,'(a,$)')'TRy this on anyway? ' read(*,'(a)')ANSWER IF (ANSWER .EQ. 'E') RETURN IF (ANSWER .NE. 'Y') GO TO 2000 END IF END IF END IF C IF (WATCH .GT. 0.5) THEN CALL SHOW (PIC(LX,LY), FMAX, SKY(ISTAR), . MX-LX+1, MY-LY+1, NCOL) WRITE (6,622) ID(ISTAR), XCEN(ISTAR), YCEN(ISTAR), . APMAG(ISTAR), NINT(FMAX) 622 FORMAT (1X, I5, 2F9.2, F9.3, 2X, 'Brightest pixel:', I7) write(*,'(a,$)')'Use this one? ' read(*,'(a)')ANSWER IF (ANSWER .EQ. 'E') RETURN IF (ANSWER .EQ. 'N') GO TO 2000 END IF C C Switch this star (at position ISTAR) with the one at position NSTAR+1 C (which is the highest star in the stack not already a PSF star). C Estimate the height of the best-fitting profile of type 1 (probably C a Gaussian) on the basis of the central pixel. C NSTAR = NSTAR+1 CALL SWAPP (ID, XCEN, YCEN, APMAG, SKY, ISTAR, NSTAR) HPSF(NSTAR) = (FMAX-SKY(NSTAR))/ . PROFIL(1, 0., 0., PAR, DFDX, DFDY, TERM, 0) SATR8D(NSTAR) = SATUR8 IF (NSTAR .LT. MAXN) GO TO 2000 C 3000 close (1) ! read *.lst end NEXP = NTAB(IEXPAND) + 2*IFRAC IF (NSTAR .LT. NEXP) THEN write(*,*) . 'There aren''t enough PSF stars for a PSF this variable.' WRITE (6,*) ' Please change something.' RETURN END IF C C If the first star is saturated, exchange it for the first C unsaturated one. C IF (SATR8D(1)) THEN DO ISTAR=2,NSTAR IF (.NOT. SATR8D(ISTAR)) GO TO 3005 END DO write(*,*)'Every single PSF star is saturated.' RETURN C 3005 CONTINUE CALL SWAPP (ID, XCEN, YCEN, APMAG, SKY, 1, ISTAR) SATR8D(ISTAR) = .TRUE. SATR8D(1) = .FALSE. END IF C CALL OVRWRT (' Chi Parameters...', 1) C C If OPT(16) (Analytic model PSF) is negative, then try all PSF types C from 1 to | OPT(16) |, inclusive. C K = NINT(OPT(16)) J = MAX0(1,K) K = IABS(K) OLD = 1.E38 IPSTYP = 0 DO 3050 I=J,K C DO L=1,NSTAR HJNK(L) = HPSF(L) XJNK(L) = XCEN(L) YJNK(L) = YCEN(L) END DO C N = NPARAM(I, FWHM, LABEL, PAR, MAXPAR) CALL FITANA (PIC, NCOL, NROW, HJNK, XJNK, YJNK, SKY, . SATR8D, NSTAR, FITRAD, WATCH, I, PAR, N, SIG) C C SIG is the root-mean-square scatter about the best-fitting analytic C function averaged over the central disk of radius FITRAD, expressed C as a fraction of the peak amplitude of the analytic model. C IF (PAR(1) .LT. 0.) GO TO 3050 C C If this latest fit is better than any previous one, store the C parameters of the stars and the model for later reference. C IF (SIG .LT. OLD) THEN IPSTYP = I WRITE (LINE,201) (PAR(L), L=1,N) 201 FORMAT (1X, 1P, 6E13.6) DO L=1,NSTAR HPSF(L) = HJNK(L) XCEN(L) = XJNK(L) YCEN(L) = YJNK(L) END DO OLD = SIG END IF 3050 CONTINUE IF (IPSTYP .LE. 0) RETURN SIG = OLD XMID = REAL(NCOL-1)/2. YMID = REAL(NROW-1)/2. xmid=0. ymid=0. do 111 i=1,nstar xmid=xmid+xcen(i) 111 ymid=ymid+ycen(i) xmid=xmid/nstar ymid=ymid/nstar NPAR = NPARAM(IPSTYP, FWHM, LABEL, PAR, MAXPAR) READ (LINE,201) (PAR(L), L=1,NPAR) C C This last bit ensures that the subsequent computations will be C performed with the parameters rounded off exactly the same as they C appear in the output .PSF file. C C======================================================================= C C At this point, we have values for the parameters of the best-fitting C analytic function. Now we may want to generate the look-up table of C corrections from the best-fitting analytic function to the actual C data. This will be generated within a square box of size NPSF x NPSF C with half-pixel spacing, centered on the centroid of the star. C C First, subtract the analytic function from all the PSF stars in the C original image. C 3300 DP = 0. SATUR8 = .FALSE. DO 3400 ISTAR=1,NSTAR RDX = AMAX1(PROFIL(IPSTYP,0.,0.5,PAR,DFDX,DFDY,TERM,0), . PROFIL(IPSTYP,0.,-0.5,PAR,DFDX,DFDY,TERM,0), . PROFIL(IPSTYP,0.5,0.,PAR,DFDX,DFDY,TERM,0), . PROFIL(IPSTYP,-0.5,0.,PAR,DFDX,DFDY,TERM,0)) IF (HPSF(ISTAR)*RDX+SKY(ISTAR) .GT. HIBAD) SATR8D(ISTAR)=.TRUE. C C The centroids and peak heights are not yet known for saturated C stars. C IF (SATR8D(ISTAR)) THEN SATUR8 = .TRUE. GO TO 3400 END IF C LX = MAX0( 1, INT(XCEN(ISTAR)-PSFRAD)-1 ) LY = MAX0( 1, INT(YCEN(ISTAR)-PSFRAD)-1 ) MX = MIN0(NCOL, INT(XCEN(ISTAR)+PSFRAD)+2 ) MY = MIN0(NROW, INT(YCEN(ISTAR)+PSFRAD)+2 ) RDX = 0.0 RDY = 0.0 DO 3350 J=LY,MY DY = REAL(J) - YCEN(ISTAR) W = DY**2 DO I=LX,MX DX = REAL(I) - XCEN(ISTAR) PIC(I,J) = PIC(I,J) - HPSF(ISTAR) * . PROFIL(IPSTYP, DX, DY, PAR, DFDX, DFDY, TERM, 0) IF (DX**2 + W .LT. RADSQ) THEN RDX = RDX + (PIC(I,J)-SKY(ISTAR))**2 RDY = RDY + 1. END IF END DO 3350 CONTINUE DP = DP + RDY ! Total number of pixels RDX = SQRT(RDX/RDY) ! Scatter inside fit radius PSF(ISTAR,1,1) = RDX/(HPSF(ISTAR)* . PROFIL(IPSTYP, 0., 0., PAR, DFDX, DFDY, TERM, 0)) 3400 CONTINUE C DP = SQRT(DP/(DP-REAL(NEXP+3*NSTAR))) ! Degrees of freedom DO ISTAR=1,NSTAR IF (.NOT. SATR8D(ISTAR)) THEN PSF(ISTAR,1,1) = DP * PSF(ISTAR,1,1) FLAG(ISTAR) = ' ' IF (SIG .GT. 0.) THEN RDX = PSF(ISTAR,1,1)/SIG IF (HPSF(ISTAR) .LE. 0) THEN WEIGHT(ISTAR) = 0. ELSE WEIGHT(ISTAR) = 1. /(1. + (RDX/2.)**2) END IF IF (RDX .GE. 3.) THEN FLAG(ISTAR) = '*' if(ifail.eq.0)ifail=1 iflag(istar)=1 ELSE IF (RDX .GE. 2.) THEN FLAG(ISTAR) = '?' END IF END IF END IF END DO K = (NSTAR-1)/5 + 1 ! Number of lines WRITE (6,67) 67 FORMAT (/' Profile errors:'/) DO I=1,K LINE = ' ' DO ISTAR=I,NSTAR,K J = 16*(ISTAR-I)/K + 1 IF (SATR8D(ISTAR)) THEN WRITE (LINE(J:J+15),68) ID(ISTAR) 68 FORMAT (1X, I5, ' saturated') ELSE WRITE (LINE(J:J+15),69) ID(ISTAR), . PSF(ISTAR,1,1), FLAG(ISTAR) 69 FORMAT (1X, I5, F7.3, 1X, A1, 1X) END IF END DO WRITE (6,*) LINE(1:79) END DO if(ifail.eq.1)return SCALE = HPSF(1) DO ISTAR=NSTAR,1,-1 HPSF(ISTAR) = HPSF(ISTAR)/HPSF(1) END DO C C Tabulate the constant part of the PSF. C MIDDLE = (NPSF+1)/2 MIDSQ = MIDDLE**2 NPSFSQ = ((NPSF-1)/2)**2 DO J=1,NPSF JDY = J-MIDDLE JDYSQ = JDY**2 DY = REAL(JDY)/2. DO I=1,NPSF IDX = I-MIDDLE K = IDX**2 + JDYSQ IF (K .LE. MIDSQ) THEN IN(I,J) = .TRUE. DX = REAL(IDX)/2. CON(I,J) = . PROFIL(IPSTYP, DX, DY, PAR, DFDX, DFDY, V, 0) IF (K .GE. NPSFSQ) EDGE(I,J) = .TRUE. ELSE IN(I,J) = .FALSE. DO K=1,NEXP PSF(I,J,K) = 0. END DO END IF END DO END DO IF (IEXPAND+IFRAC .LT. 0) THEN PSFMAG = 0.0D0 DO J=1,NPSF JDY = J-MIDDLE JDYSQ = JDY**2 DO I=1,NPSF IF (IN(I,J)) THEN IDX = I-MIDDLE IF (IDX**2+JDYSQ .LE. IRADSQ) PSFMAG = PSFMAG + . CON(I,J) END IF END DO END DO PSFMAG = 25.D0-2.5D0*DLOG10(SCALE*PSFMAG/4.D0) GO TO 2900 END IF C C======================================================================= C C Now compute look-up tables. C 3500 CONTINUE C C MIDDLE is the center of the look-up table, which will correspond to C the centroid of the analytic PSF. MIDSQ is the square of the C radius of the PSF table, plus a pixel's worth of slack just to C make sure you'll always have a 4x4 array suitable for interpolation. C DO K=1,NEXP DO J=1,NPSF DO I=1,NPSF PSF(I,J,K) = 0.0 END DO END DO END DO C NPASS = AMIN0(NSTAR-1,NPASS) DO 4950 JY=1,NPSF IF (WATCH .GT. -1.5) THEN WRITE (LINE,*) ' Computed', JY, ' rows of', NPSF, . ' in the PSF.' CALL OVRWRT (LINE(1:60), 2) END IF RDY = REAL(JY-MIDDLE)/2. C DO 4940 IX=1,NPSF C C Don't waste time with pixels outside a radius = (MIDDLE-1). C IF (.NOT. IN(IX,JY)) GO TO 4940 RDX = REAL(IX-MIDDLE)/2. DO 4935 ITER=0,NPASS C C Initialize accumulators. C TERM(1) = 1. SUMSQ = 0. SUMN = 0. DO L=1,NEXP V(L) = 0. DO K=1,NEXP C(K,L) = 0. END DO END DO C C Now, for this point in the set of lookup table(s) [(IX,JY), where C (MIDDLE,MIDDLE) corresponds to the centroid of the PSF], consult C each of the PSF stars to determine the value(s) to put into the C table(s). C DO 4900 ISTAR = 1,NSTAR IF (SATR8D(ISTAR)) GO TO 4900 C C What pixels in the original image constitute a 4x4 box surrounding C this (IX,JY) in the PSF, as referred to this star's centroid? C DX = XCEN(ISTAR) + RDX DY = YCEN(ISTAR) + RDY I = INT(DX) J = INT(DY) C C The 4x4 box is given by I-1 <= x <= I+2, J-1 <= y <= J+2. C IF ((I .LT. 2) .OR. (J .LT. 2) .OR. . (I+2 .GT. NCOL) .OR. . (J+2 .GT. NROW)) GO TO 4900 ! Next star C DO L=J-1,J+2 DO K=I-1,I+2 IF ((PIC(K,L) .GT. HIBAD) .OR. . (PIC(K,L) .LT. -HIBAD)) GO TO 4900 END DO END DO C DX = DX - I DY = DY - J C C The point which corresponds PRECISELY to the offset (RDX,RDY) C from the star's centroid, lies a distance (DX,DY) from pixel C (I,J). C C Use bicubic interpolation to evaluate the residual PSF amplitude C at this point. Scale the residual up to match the first PSF star. C DP = BICUBC(PIC(I-1,J-1), NCOL, DX, DY, DFDX, DFDY) . - SKY(ISTAR) DP = DP/HPSF(ISTAR) SUMSQ = SUMSQ + ABS(DP) SUMN = SUMN + 1. IF (IEXPAND .GE. 1) THEN TERM(2) = (XCEN(ISTAR)-1.)/XMID-1. TERM(3) = (YCEN(ISTAR)-1.)/YMID-1. IF (IEXPAND .GE. 2) THEN TERM(4) = 1.5*TERM(2)**2-0.5 TERM(5) = TERM(2)*TERM(3) TERM(6) = 1.5*TERM(3)**2-0.5 IF (IEXPAND .GE. 3) THEN TERM(7) = TERM(2)*(5.*TERM(4)-2.)/3. TERM(8) = TERM(4)*TERM(3) TERM(9) = TERM(2)*TERM(6) TERM(10) = TERM(3)*(5.*TERM(6)-2.)/3. END IF END IF END IF C C IF (IFRAC .GE. 1) THEN C C INSERT CODE HERE C IF (ITER .GT. 0) THEN OLD = 0. DO K=1,NEXP OLD = OLD + PSF(IX,JY,K) * TERM(K) END DO C IF (ITER .LE. MAX0(3,NPASS/2)) THEN W = HPSF(ISTAR)/(1. + (ABS(DP - OLD)/SIG)) ELSE W = HPSF(ISTAR)/(1. + ((DP - OLD)/SIG)**2) END IF ELSE W = HPSF(ISTAR) END IF C W = W*WEIGHT(ISTAR) DO K=1,NEXP WT = W*TERM(K) V(K) = V(K) + WT*DP DO L=1,NEXP C(K,L) = C(K,L) + WT*TERM(L) END DO END DO 4900 CONTINUE C IF (SUMN .LT. NEXP) THEN write(*,*)'Not enough PSF stars. Please start over.' close(3) RETURN END IF CALL INVERS (C, MAXEXP, NEXP, ISTAT) CALL VMUL (C, MAXEXP, NEXP, V, A) DO K=1,NEXP PSF(IX,JY,K) = A(K) END DO IF (SUMN .LE. NEXP) GO TO 4940 SIG = 1.2533*SUMSQ/SQRT(SUMN*(SUMN - NEXP)) 4935 CONTINUE 4940 CONTINUE 4950 CONTINUE C C Make the average value outside the PSF radius identically zero. C DX = 0. K = 0 DO J=1,NPSF DO I=1,NPSF IF (EDGE(I,J)) THEN DX = DX+SCALE*CON(I,J)+PSF(I,J,1) K = K+1 END IF END DO END DO DX = DX/REAL(K) DO J=1,NPSF DO I=1,NPSF IF (IN(I,J)) PSF(I,J,1) = PSF(I,J,1)-DX END DO END DO IF (NEXP .LT. 2) GO TO 2800 C C At this point, we must be sure that any higher order terms C in the PSF [i.e. those that go as powers of (XCEN-XMID) and C (YCEN-YMID)] contain zero volume, so that the total volume C of the PSF is independent of position. This will be done C by looking at the total flux contained in the look-up C tables of index 2 and higher. The analytic function C will be scaled to the same total flux and subtracted from C each of the higher lookup tables and added into the C lookup table of index 1. I scale the analytic function C instead of simply transferring the net flux from the C higher tables to table 1, because it is poor fits of C the analytic profile to the image data which has caused C the net flux in the lookup tables to depend upon position. C C Compute the brightness weighted average value of each of C the terms in the polynomial expansion. C DO K=1,NEXP TERM(K) = 0. END DO C DO ISTAR=NSTAR,1,-1 DX = (XCEN(ISTAR)-1.)/XMID-1. DY = (YCEN(ISTAR)-1.)/YMID-1. W = WEIGHT(ISTAR)*HPSF(ISTAR) TERM(1) = TERM(1) + W TERM(2) = TERM(2) + W*DX TERM(3) = TERM(3) + W*DY IF (IEXPAND .GE. 2) THEN DX2 = DX**2 DY2 = DY**2 TERM(4) = TERM(4) + W*(1.5*DX2-0.5) TERM(5) = TERM(5) + W*(DX*DY) TERM(6) = TERM(6) + W*(1.5*DY2-0.5) IF (IEXPAND .GE. 3) THEN TERM(7) = TERM(7) + DX*(5.*DX2-3.)/2. TERM(8) = TERM(8) + DY*(1.5*DX2-0.5) TERM(9) = TERM(9) + DX*(1.5*DY2-0.5) TERM(10) = TERM(10) + DY*(5.*DY2-3.)/2. END IF END IF C C IF (IFRAC .GE. 1) THEN C C INSERT CODE HERE C END DO C DO K=NEXP,1,-1 TERM(K) = TERM(K)/TERM(1) END DO C C Okey doke. Now if there is any net volume contained in any of the C higher order (variable) terms of the PSF, we will remove it in two C components. First of all, there is the possibility that there were some C errors in the sky estimates which was correlated with position. C Then there is the possibility that there was some systematic error C in estimating HPSF across the frame. We correct these by removing C any net volume in the higer-order look-up tables and by removing C any correlation with the analytic profile. We want to make sure C we do this only for the area within one PSF radius. C C DO K=2,NEXP DO JJ=0,1 DO II=0,1 DO J=1,2 VV(J) = 0.D0 DO I=1,2 CC(I,J) = 0.D0 END DO END DO DO J=1+JJ,NPSF-JJ,2 JDY = J-MIDDLE JDYSQ = JDY**2 DO I=1+II,NPSF-II,2 IF (IN(I,J)) THEN IDX = I-MIDDLE IR = IDX**2+JDYSQ IF (IR .LE. IRADSQ) THEN VV(1) = VV(1)+PSF(I,J,K) cxyz VV(2) = VV(2)+PSF(I,J,K)*CON(I,J) CC(1,1) = CC(1,1)+1. CC(2,1) = CC(2,1)+CON(I,J) cxyz CC(2,2) = CC(2,2)+CON(I,J)**2 END IF IF (IR .GE. JRADSQ) THEN VV(2) = VV(2)+PSF(I,J,K) CC(1,2) = CC(1,2)+1. CC(2,2) = CC(2,2)+CON(I,J) END IF END IF END DO END DO cxyz CC(1,2) = CC(2,1) CALL DINVRS (CC, 2, 2, ISTAT) CALL DVMUL (CC, 2, 2, VV, AA) SCALE = SCALE + TERM(K)*AA(2) DX = TERM(K)*AA(1) DO J=1+JJ,NPSF-JJ,2 JDY = J-MIDDLE JDYSQ = JDY**2 DO I=1+II,NPSF-II,2 IF (IN(I,J)) THEN IDX = I-MIDDLE PSF(I,J,K) = PSF(I,J,K) - . AA(1)-CON(I,J)*AA(2) PSF(I,J,1) = PSF(I,J,1) + DX END IF END DO END DO END DO END DO END DO C 2799 continue C C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! C 2800 CONTINUE PSFMAG = 0.0D0 DO J=1,NPSF JDY = J-MIDDLE JDYSQ = JDY**2 DO I=1,NPSF IF (IN(I,J)) THEN IDX = I-MIDDLE IF (IDX**2+JDYSQ .LE. IRADSQ) THEN PSFMAG = PSFMAG + SCALE*CON(I,J) + DBLE(PSF(I,J,1)) END IF END IF END DO END DO if(PSFMAG.lt.0)then inquire(file='psf.mag',exist=logi) if(.not.logi)then PSFMAG=12.5 else open(55,file='psf.mag',status='old') read(55,*)PSFMAG PSFMAG=PSFMAG+0.5 close(55) endif else PSFMAG = 25.D0-2.5D0*DLOG10(PSFMAG/4.D0) open(55,file='psf.mag',status='unknown') write(55,*)PSFMAG close(55) endif C C If there are saturated stars, we must now guesstimate their positions C and peak brightnesses. Determine these from single-profile fits a la C the PEAK routine, with a fitting radius equal to half the PSF radius. C Subtract the analytic part of their profile from the image, and C go back and redetermine the lookup tables. C IF (SATUR8 .AND. (OPT(18) .GT. 0.5)) THEN PERR = 0.01*OPT(19) PKERR = 0.01*OPT(20)/(PAR(1)*PAR(2)) ! **2 DO ISTAR=1,NSTAR IF (SATR8D(ISTAR)) THEN LX = MAX0( 1, INT(XCEN(ISTAR)-PSFRAD)-1 ) LY = MAX0( 1, INT(YCEN(ISTAR)-PSFRAD)-1 ) MX = MIN0(NCOL, INT(XCEN(ISTAR)+PSFRAD)+2 ) MY = MIN0(NROW, INT(YCEN(ISTAR)+PSFRAD)+2 ) K = MX-LX+1 L = MY-LY+1 DX = (XCEN(ISTAR)-1.)/XMID-1. DY = (YCEN(ISTAR)-1.)/YMID-1. HPSF(ISTAR) = 3.*HPSF(1) ! Lousy starting guess XCEN(ISTAR) = XCEN(ISTAR)-LX+1. YCEN(ISTAR) = YCEN(ISTAR)-LY+1. CALL PKFIT (PIC(LX,LY), K, L, NCOL, XCEN(ISTAR), . YCEN(ISTAR), HPSF(ISTAR), SKY(ISTAR), 0.5*PSFRAD, . LOBAD, HIBAD, SCALE, IPSTYP, PAR, MAXPAR, NPAR, . PSF, MAXPSF, MAXEXP, NPSF, NEXP, IFRAC, DX, . DY, ERRMAG, CHI, SHARP, N, NCOL, NROW) XCEN(ISTAR) = XCEN(ISTAR)+LX-1. YCEN(ISTAR) = YCEN(ISTAR)+LY-1. WRITE (6,666) ID(ISTAR), XCEN(ISTAR), YCEN(ISTAR), . PSFMAG-2.5*ALOG10(HPSF(ISTAR)), . AMIN1(2.0, 1.086*ERRMAG/HPSF(ISTAR)), . REAL(N), CHI, SHARP 666 FORMAT (1X, I5, 4F9.3, F9.0, F9.2, F9.3) LX = MAX0( 1, INT(XCEN(ISTAR)-PSFRAD)-1 ) LY = MAX0( 1, INT(YCEN(ISTAR)-PSFRAD)-1 ) MX = MIN0(NCOL, INT(XCEN(ISTAR)+PSFRAD)+2 ) MY = MIN0(NROW, INT(YCEN(ISTAR)+PSFRAD)+2 ) RDX = SCALE*HPSF(ISTAR) DO J=LY,MY DY = REAL(J) - YCEN(ISTAR) W = DY**2 DO I=LX,MX IF (PIC(I,J) .LE. HIBAD) THEN DX = REAL(I) - XCEN(ISTAR) RDY = RDX * PROFIL(IPSTYP, . DX, DY, PAR, DFDX, DFDY, TERM, 0) IF (RDY+SKY(ISTAR) .LE. HIBAD) THEN PIC(I,J) = PIC(I,J) - RDY ELSE PIC(I,J) = RDY+SKY(ISTAR) C C Ensuring it is greater than HIBAD, and hence will be ignored. C END IF END IF END DO END DO SATR8D(ISTAR) = .FALSE. WEIGHT(ISTAR) = 0.5 END IF END DO SATUR8 = .FALSE. GO TO 3500 END IF C C--------------------------------------------------------------------- C 2900 CONTINUE C C Write PSF file header. C WRITE (3,202) LABEL, NPSF, NPAR, NTAB(IEXPAND), 5*IFRAC, . PSFMAG, SCALE, XMID, YMID WRITE (*,202) LABEL, NPSF, NPAR, NTAB(IEXPAND), 5*IFRAC, . PSFMAG, SCALE, XMID, YMID 202 FORMAT (1X, A8, 4I5, F9.3, F15.3, 2F9.1) WRITE (3,201) (PAR(L), L=1,NPAR) C IF (NEXP .GT. 0) THEN DO K=1,NEXP WRITE (3,204) ((PSF(I,J,K), I=1,NPSF), J=1,NPSF) 204 FORMAT (1X, 1P, 6E13.6) END DO END IF close(3) C C Now generate the file of neighbors. C C First, write out the PSF stars. C CALL WRHEAD (4, 3, NCOL, NROW, 7, LOBAD, HIBAD, THRESH, AP1, . PHPADU, READNS, RADIUS) DO I=1,NSTAR WRITE (4,280) ID(I), XCEN(I), YCEN(I), APMAG(I), SKY(I) 280 FORMAT (1X, I5, 3F9.3, f11.3) END DO IF (NSTAR .GE. NTOT) GO TO 9000 C C Now look for all stars within a radius equal to C 1.5*PSFRAD+2.*FITRAD+1 of any of the PSF stars. C RSQ = (1.5*PSFRAD+2.*FITRAD+1.)**2 L = NSTAR+1 C C L will point at the first irrelevant star. C DO I=1,NSTAR J = L 5020 IF ((XCEN(J)-XCEN(I))**2 + (YCEN(J)-YCEN(I))**2 .LE. RSQ) THEN CALL SWAPP (ID, XCEN, YCEN, APMAG, SKY, J, L) WRITE (4,280) ID(L), XCEN(L), YCEN(L), APMAG(L), SKY(L) L = L+1 END IF IF (J .LT. NTOT) THEN J = J+1 GO TO 5020 END IF END DO C C Now all stars in the stack between positions NSTAR+1 and L-1, C inclusive, lie within the specified radius of the PSF stars and have C been written out. Finally, look for any stars within a radius C equal to 2*FITRAD+1 of any of THESE stars, and of each other. C IF (L .GT. NTOT) GO TO 9000 RSQ = (2.*FITRAD+1)**2 I = NSTAR+1 5110 J = L 5120 IF ((XCEN(J)-XCEN(I))**2 + (YCEN(J)-YCEN(I))**2 .LE. RSQ) THEN CALL SWAPP (ID, XCEN, YCEN, APMAG, SKY, J, L) WRITE (4,280) ID(L), XCEN(L), YCEN(L), APMAG(L), SKY(L) L = L+1 END IF IF (J .LT. NTOT) THEN J = J+1 GO TO 5120 END IF C I = I+1 IF (I .LT. L) GO TO 5110 9000 close(4) RETURN END! C C####################################################################### C SUBROUTINE FITANA (PIC, NCOL, NROW, H, XCEN, YCEN, SKY, . SATR8D, NSTAR, . FITRAD, WATCH, IPSTYP, PAR, NPAR, CHI) C C This subroutine fits the ANALYTIC profile to the selected stars. C C OFFICIAL DAO VERSION: 1991 May 10 C IMPLICIT NONE INTEGER MAXPAR, MAXBOX, MAXN, NCOL, NROW PARAMETER (MAXPAR=6, MAXBOX=69, MAXN=200) CHARACTER*80 LINE REAL PIC(NCOL,NROW) REAL H(*), XCEN(*), YCEN(*), SKY(*), CLAMP(MAXPAR), OLD(MAXPAR) REAL C(MAXPAR,MAXPAR), V(MAXPAR), PAR(MAXPAR), T(MAXPAR), . Z(MAXPAR) REAL XCLAMP(MAXN), YCLAMP(MAXN), XOLD(MAXN), YOLD(MAXN) REAL ABS, PROFIL LOGICAL SATR8D(*) C REAL DHN, DHD, DXN, DXD, DYN, DYD, DX, DY, DYSQ REAL PEAK, DP, PROD, DHDXC, DHDYC, P, WT REAL RSQ, SUMWT, OLDCHI, CHI, WATCH, FITRAD INTEGER I, J, K, L, LX, LY, MX, MY, ISTAT INTEGER ISTAR, MPAR, NITER, NPAR, IPSTYP, NSTAR LOGICAL FULL C C----------------------------------------------------------------------- C FULL = .FALSE. DO I=1,NPAR Z(I) = 0. OLD(I) = 0. CLAMP(I) = 0.5 END DO CLAMP(1) = 2. CLAMP(2) = 2. NITER = 0 OLDCHI = 0. DO I=1,NSTAR XOLD(I) = 0. YOLD(I) = 0. XCLAMP(I) = 1. YCLAMP(I) = 1. END DO C MPAR = 2 C C----------------------------------------------------------------------- C C SECTION 1 C C Now we will fit an integrated analytic function to the central part C of the stellar profile. For each star we will solve for three C parameters: (1) H, the central height of the model profile (above C sky); (2) XCEN, C the centroid of the star in x; and (3) YCEN, C likewise for y. In addition, from ALL STARS CONSIDERED TOGETHER C we will determine any other parameters, PAR(i), required to describe C the profile. We will use a circle of radius FITRAD centered on the C position of each PSF star. NOTE THAT we do not fit the data to an C actual analytic profile, but rather the function is numerically C integrated over the area of each pixel, and the observed data are fit C to these integrals. C RSQ = FITRAD**2 1000 NITER = NITER+1 IF (NITER .GT. 600) THEN write(*,*)'Failed to converge.' PAR(1) = -1. RETURN END IF C C Initialize the big accumulators. C DO I=1,NPAR V(I)=0.0 ! Zero the vector of residuals DO J=1,NPAR C(J,I)=0.0 ! Zero the normal matrix END DO END DO CHI = 0.0 SUMWT = 0.0 C C Using the analytic model PSF defined by the current set of parameters, C compute corrections to the brightnesses and centroids of all the PSF C stars. MEANWHILE, accumulate the corrections to the model parameters. C DO 1400 ISTAR=1,NSTAR IF (SATR8D(ISTAR)) GO TO 1400 LX = INT(XCEN(ISTAR)-FITRAD) + 1 LY = INT(YCEN(ISTAR)-FITRAD) + 1 MX = INT(XCEN(ISTAR)+FITRAD) MY = INT(YCEN(ISTAR)+FITRAD) C DHN = 0.0D0 DHD = 0.0D0 DXN = 0.0D0 DXD = 0.0D0 DYN = 0.0D0 DYD = 0.0D0 DO J=LY,MY DY = REAL(J) - YCEN(ISTAR) DYSQ = DY**2 DO I=LX,MX DX = REAL(I) - XCEN(ISTAR) WT = (DX**2+DYSQ)/RSQ IF (WT .LT. 1.) THEN P = PROFIL(IPSTYP, DX, DY, PAR, DHDXC, DHDYC, T, 0) DP = PIC(I,J) - H(ISTAR)*P - SKY(ISTAR) c if (istar .le. 2) type 6661, i, j, pic(i,j), c . h(istar)*p, sky(istar), dp c6661 format (2i6, 4f12.3) DHDXC = H(ISTAR)*DHDXC DHDYC = H(ISTAR)*DHDYC WT = 5./(5.+WT/(1.-WT)) PROD = WT*P DHN = DHN + PROD*DP DHD = DHD + PROD*P PROD = WT*DHDXC DXN = DXN + PROD*DP DXD = DXD + PROD*DHDXC PROD = WT*DHDYC DYN = DYN + PROD*DP DYD = DYD + PROD*DHDYC END IF END DO END DO C c if (istar .le. 2) type 6662, h(istar), dhn, dhd, c . h(istar)+dhn/dhd c6662 format (4f12.3) c if (istar .eq. 2) type * c if (istar .eq. 2) accept * H(ISTAR) = H(ISTAR) + DHN/DHD C DXN = DXN/DXD IF (XOLD(ISTAR)*DXN .LT. 0.) XCLAMP(ISTAR) = 0.5*XCLAMP(ISTAR) XOLD(ISTAR) = DXN XCEN(ISTAR) = XCEN(ISTAR)+DXN/(1.+ABS(DXN)/XCLAMP(ISTAR)) C DYN = DYN/DYD IF (YOLD(ISTAR)*DYN .LT. 0.) YCLAMP(ISTAR) = 0.5*YCLAMP(ISTAR) YOLD(ISTAR) = DYN YCEN(ISTAR) = YCEN(ISTAR)+DYN/(1.+ABS(DYN)/YCLAMP(ISTAR)) C PEAK = H(ISTAR) * PROFIL(IPSTYP,0.,0., PAR, DHDXC, DHDYC, T, 0) DO J=LY,MY DY = REAL(J)-YCEN(ISTAR) DYSQ = DY**2 DO I=LX,MX DX = REAL(I)-XCEN(ISTAR) WT = (DX**2+DYSQ)/RSQ IF (WT .LT. 1.) THEN P = PROFIL(IPSTYP, DX, DY, PAR, DHDXC, DHDYC, T, 1) DP = PIC(I,J) - H(ISTAR)*P - SKY(ISTAR) DO K=1,MPAR T(K) = H(ISTAR)*T(K) END DO CHI = CHI + ABS(DP/PEAK) SUMWT = SUMWT + 1. WT = 5./(5.+WT/(1.-WT)) IF (NITER .GE. 4) WT = WT/(1.+ABS(20.*DP/PEAK)) DO K=1,MPAR V(K) = V(K) + WT*DP*T(K) DO L=1,MPAR C(L,K) = C(L,K) + WT*T(L)*T(K) END DO END DO END IF END DO END DO 1400 CONTINUE C C Correct the fitting parameters. C CALL INVERS (C, MAXPAR, MPAR, ISTAT) CALL VMUL (C, MAXPAR, MPAR, V, Z) DO I=1,MPAR IF (Z(I)*OLD(I) .LT. 0.) THEN CLAMP(I) = 0.5*CLAMP(I) ELSE CLAMP(I) = 1.1*CLAMP(I) END IF OLD(I) = Z(I) Z(I) = CLAMP(I)*Z(I) END DO Z(1) = AMAX1(-0.1*PAR(1), AMIN1(0.1*PAR(1), Z(1))) Z(2) = AMAX1(-0.1*PAR(2), AMIN1(0.1*PAR(2), Z(2))) C Z(3) = Z(3)/(1.+ABS(Z(3))/(AMIN1(0.1,1.-ABS(PAR(3))))) DO I=1,MPAR PAR(I) = PAR(I)+Z(I) END DO C SUMWT = SUMWT - REAL(MPAR + 3*NSTAR) IF (SUMWT .GT. 0) THEN CHI = 1.2533*CHI/SUMWT ELSE CHI = 9.9999 END IF C WRITE (LINE,661) CHI, (PAR(I), I=1,MPAR) 661 FORMAT (1X, F7.4, 6F10.5) I=10*MPAR+13 IF (WATCH .GT. -1.5) CALL OVRWRT (LINE(1:I), 2) IF (MPAR .EQ. NPAR) THEN IF (ABS(OLDCHI/CHI-1.) .LT. 1.E-5) THEN CALL OVRWRT (LINE(1:I), 3) RETURN END IF ELSE IF (ABS(OLDCHI/CHI-1.) .LT. 1.E-3) THEN MPAR = MPAR+1 OLDCHI = 0. GO TO 1000 END IF END IF C OLDCHI = CHI GO TO 1000 END! C C####################################################################### C C Exchange two stars in the star list. C C OFFICIAL DAO VERSION: 1991 May 10 C SUBROUTINE SWAPP (ID, X, Y, A, S, I, J) IMPLICIT NONE REAL X(*), Y(*), A(*), S(*) INTEGER ID(*) C REAL HOLD INTEGER I, J, IHOLD C IHOLD = ID(I) ID(I) = ID(J) ID(J) = IHOLD HOLD = X(I) X(I) = X(J) X(J) = HOLD HOLD = Y(I) Y(I) = Y(J) Y(J) = HOLD HOLD = A(I) A(I) = A(J) A(J) = HOLD HOLD = S(I) S(I) = S(J) S(J) = HOLD RETURN END! C C####################################################################### C SUBROUTINE WRHEAD (LUN, NL, NCOL, NROW, ITEMS, LOBAD, HIBAD, . THRESH, AP1, PHPADU, READNS, FRAD) C C======================================================================= C C Subroutine to write a standard header into an output sequential C data file on the disk. Same as RDHEAD, except that all the output C arguments are now input arguments. ITEMS tells how many of the C individual arguments are to be written into the header. C C Latest DAOPHOT II version: 1991 January 15 C======================================================================= C IMPLICIT NONE CHARACTER*8 HEAD(7) REAL A(7) C REAL LOBAD, HIBAD, THRESH, AP1, PHPADU, READNS, FRAD INTEGER LUN, NL, NCOL, NROW, ITEMS, I DATA HEAD /' LOWBAD', ' HIGHBAD', ' THRESH', ' AP1', . ' PH/ADU', ' RNOISE', ' FRAD'/ C C----------------------------------------------------------------------- C A(1)=LOBAD A(2)=HIBAD A(3)=THRESH A(4)=AP1 A(5)=PHPADU A(6)=READNS A(7)=FRAD WRITE (LUN,900) (HEAD(I), I=1,ITEMS) 900 FORMAT (' NL NX NY', 8A8) IF (A(2) .LT. 99999.9) THEN WRITE (LUN,901) NL, NCOL, NROW, (A(I), I=1,ITEMS) 901 FORMAT (1X, I2, 2I5, 2F8.1, 5F8.2) ELSE WRITE (LUN,902) NL, NCOL, NROW, (A(I), I=1,ITEMS) 902 FORMAT (1x, I2, 2I5, F8.0, F9.0, 5F8.2) END IF WRITE (LUN,901) ! Write a blank line RETURN ! Normal return C END! C C####################################################################### C SUBROUTINE RDSTAR (LU, NL, ID, X, Y, AMAG, SKY) C C======================================================================= C C Read in an ID number, x and y coordinates, magnitude, and sky value C for the next star in the input file, whatever the file type (i.e., C NL = 1, 2, or 3). If an end of file is encountered, set the star ID C to a negative number and return; if a blank line is encountered, the C star ID will be zero. LU is the logical unit number to be used; the C other arguments are obvious. C C======================================================================= C IMPLICIT NONE CHARACTER*133 LINE REAL X, Y, AMAG, SKY, DUMMY INTEGER LU, N, ISTAT, ID, NL C 900 CALL RDCHAR (LU, LINE, N, ISTAT) IF (ISTAT .GT. 0) THEN ID = -1 ! ID = -1 means END OF FILE RETURN ELSE IF (ISTAT .LT. 0) THEN write(*,*)'Unable to read line from input file' END IF C IF (N .LE. 1) THEN ID = 0 ! Blank line RETURN END IF C N = N+1 LINE(N:N) = '/' IF (NL .EQ. 1) THEN READ (LINE(1:N),*,ERR=1000) ID, X, Y, AMAG, DUMMY, SKY ELSE IF (NL .EQ. 2) THEN READ (LINE(1:N),*,ERR=1000,END=2000) ID, X, Y, AMAG IF (ID .NE. 0) READ(LU,321,ERR=1000) SKY 321 FORMAT (4X, F9.3) ELSE IF (NL .EQ. 3) THEN READ (LINE(1:N),*,ERR=1000) ID, X, Y, AMAG, SKY END IF RETURN ! Normal return 1000 write(*,*)'WARNING: Corrupt star data encountered in input:' GO TO 900 2000 ID=-1 RETURN END! C C####################################################################### C SUBROUTINE RDCHAR (LUN, LINE, N, ISTAT) IMPLICIT NONE CHARACTER*(*) LINE INTEGER LUN, N, ISTAT, I, J C READ (LUN,1,END=8000,ERR=9000) LINE 1 FORMAT (A) N = 0 DO I=1,LEN(LINE) J = ICHAR(LINE(I:I)) IF ((J .GE. 33) .AND. (J .LE. 126)) N = I END DO ISTAT = 0 RETURN 8000 LINE = ' ' N = 0 ISTAT = 1 ! ISTAT = 1 means END OF FILE RETURN 9000 LINE = ' ' N = 0 ISTAT = -1 ! ISTAT = -1 means ERROR RETURN END! C C This file contains subroutines that are not I/O related, but rather C involve arithmetic or character operations that my be somewhat C machine-specific-- it may be necessary or desirable to make some C changes to optimize the code to run on your computer. C C OFFICIAL DAO VERSION: 1991 April 18 C C*********************************************************************** C C Current contents: C C INVERS inverts a square matrix. C C VMUL multiplies a square matrix (on the left) by a column vector C (on the right) yielding a new column vector. C C QUICK does a quicksort on a vector of data. C C SHOW produces what is effectively an eleven-level gray scale plot C of a rectangular data array on the terminal. C C*********************************************************************** C SUBROUTINE INVERS (A, MAX, N, IFLAG) C C Although it seems counter-intuitive, the tests that I have run C so far suggest that the 180 x 180 matrices that NSTAR needs can C be inverted with sufficient accuracy if the elements are REAL C rather than DOUBLE PRECISION C C Arguments C C A (INPUT/OUTPUT) is a square matrix of dimension N. The inverse C of the input matrix A is returned in A. C C MAX (INPUT) is the size assigned to the matrix A in the calling C routine. It's needed for the dimension statement below. C C IFLAG (OUTPUT) is an error flag. IFLAG = 1 if the matrix could not C be inverted; IFLAG = 0 if it could. C IMPLICIT NONE INTEGER MAX REAL A(MAX,MAX) C INTEGER N, IFLAG, I, J, K C C----------------------------------------------------------------------- C IFLAG=0 I=1 300 IF(A(I,I).EQ.0.0E0)GO TO 9100 A(I,I)=1.0E0/A(I,I) J=1 301 IF(J.EQ.I)GO TO 304 A(J,I)=-A(J,I)*A(I,I) K=1 302 IF(K.EQ.I)GO TO 303 A(J,K)=A(J,K)+A(J,I)*A(I,K) 303 IF(K.EQ.N)GO TO 304 K=K+1 GO TO 302 304 IF(J.EQ.N)GO TO 305 J=J+1 GO TO 301 305 K=1 306 IF(K.EQ.I)GO TO 307 A(I,K)=A(I,K)*A(I,I) 307 IF(K.EQ.N)GO TO 308 K=K+1 GO TO 306 308 IF(I.EQ.N)RETURN ! Normal return I=I+1 GO TO 300 C C----------------------------------------------------------------------- C C Error: zero on the diagonal. C 9100 IFLAG=1 RETURN C END! C C####################################################################### C SUBROUTINE VMUL (A, MAX, N, V, X) C C Multiply a matrix by a vector: C C A * V = X C C Arguments C C A(column,row) (INPUT) is a square matrix of dimension N. C C MAX (INPUT) is the size assigned to the array in the C calling routine. C C V(row) (INPUT) is a column vector of dimension N. C C X(row) (OUTPUT) is a column vector of dimension N. C IMPLICIT NONE INTEGER MAX REAL A(MAX,MAX), V(MAX) REAL X(MAX) C DOUBLE PRECISION DBLE REAL SNGL C DOUBLE PRECISION SUM INTEGER N, I, J C C----------------------------------------------------------------------- C I=1 200 SUM=0.0D0 J=1 201 SUM=SUM+DBLE(A(J,I))*DBLE(V(J)) IF (J .EQ. N) GO TO 203 J=J+1 GO TO 201 203 X(I)=SNGL(SUM) IF (I .EQ. N) RETURN ! Normal return I=I+1 GO TO 200 END! C C####################################################################### C SUBROUTINE DINVRS (A, MAX, N, IFLAG) IMPLICIT NONE INTEGER MAX DOUBLE PRECISION A(MAX,MAX) C INTEGER I, J, K, N, IFLAG C C----------------------------------------------------------------------- C IFLAG=0 I=1 300 IF(A(I,I).EQ.0.0D0)GO TO 9100 A(I,I)=1.0D0/A(I,I) J=1 301 IF(J.EQ.I)GO TO 304 A(J,I)=-A(J,I)*A(I,I) K=1 302 IF(K.EQ.I)GO TO 303 A(J,K)=A(J,K)+A(J,I)*A(I,K) 303 IF(K.EQ.N)GO TO 304 K=K+1 GO TO 302 304 IF(J.EQ.N)GO TO 305 J=J+1 GO TO 301 305 K=1 306 IF(K.EQ.I)GO TO 307 A(I,K)=A(I,K)*A(I,I) 307 IF(K.EQ.N)GO TO 308 K=K+1 GO TO 306 308 IF(I.EQ.N)RETURN ! Normal return I=I+1 GO TO 300 C C----------------------------------------------------------------------- C C Error: zero on the diagonal. C 9100 IFLAG=1 RETURN C END! C C####################################################################### C SUBROUTINE DVMUL (A, MAX, N, V, X) IMPLICIT NONE INTEGER MAX DOUBLE PRECISION A(MAX,MAX), V(MAX), X(MAX) C DOUBLE PRECISION SUM INTEGER I, J, N C C----------------------------------------------------------------------- C I=1 200 SUM=0.0D0 J=1 201 SUM=SUM+A(J,I)*V(J) IF (J .EQ. N) GO TO 203 J=J+1 GO TO 201 203 X(I)=SUM IF (I .EQ. N) RETURN ! Normal return I=I+1 GO TO 200 END! C C####################################################################### C SUBROUTINE SHOW (F, FMAX, FZERO, NX, NY, NCOL) C C======================================================================= C C A simple subroutine to use alphanumeric characters to produce a 2-D C gray-scale plot on an alphanumeric terminal. C C======================================================================= C IMPLICIT NONE INTEGER NCHAR, MAXPLT, NCOL PARAMETER (NCHAR=11, MAXPLT=36) C C Parameters C C NCHAR is the number of discrete gray levels we wish to produce. C C MAXPLT is the widest plot that can be produced on the terminal C screen. Since two characters will be typed out per pixel (to C make the overall plot more nearly square) MAXPLT should be C equal to (N-2)/2 where N is the number of character positions C per line on the screen; the extra two characters will be used C for vertical bars ('|') to delimit the picture. Arrays C which are more than MAXPLT pixels on a side will be C rebinned before display. C REAL F(NCOL,*), FF(MAXPLT) C REAL SQRT INTEGER MIN0 C CHARACTER BLANK*80, DASH*2 CHARACTER*2 CHAR(NCHAR) REAL SUM, PIXELS, FMAX, FZERO, S INTEGER I, ISTEP, NX, MX, IX, IY, KX, JX, JY, NY, LOW DATA BLANK /' '/, DASH/'--'/ DATA CHAR / ' ', '- ', '--', '::', '==', 'll', . 'II', '%%', '00', 'HH', '##' / C C----------------------------------------------------------------------- C ISTEP=((NX-1)/MAXPLT)+1 C C ISTEP is the number of pixels in each row of the input array which C will have to be averaged for each pixel of the display. C MX=(NX+ISTEP-1)/ISTEP C C MX is the number of pixels per row which will be produced on the C output display. C LOW=MAXPLT-MX+1 S=SQRT(AMAX1(FLOAT(NCHAR), FMAX-FZERO)) WRITE (6,610) BLANK(1:LOW), (DASH, I=1,MX), '+ ' 610 FORMAT (A, '+', 80A2) C DO 1010 IY=1,NY,ISTEP KX=0 DO 1007 IX=1,NX,ISTEP KX=KX+1 PIXELS=0.0 SUM=0.0 DO 1005 JY=IY,MIN0(NY,IY+ISTEP-1) DO 1003 JX=IX,MIN0(NX,IX+ISTEP-1) PIXELS=PIXELS+1.0 1003 SUM=SUM+F(JX,JY) 1005 CONTINUE 1007 FF(KX)=SUM/PIXELS 1010 WRITE (6,602) BLANK(1:LOW), '|', (CHAR(MIN0( . NCHAR, . IFIX( NCHAR*SQRT(AMAX1(0., FF(IX)-FZERO) )/S )+1)), . IX=1,MX), '|' 602 FORMAT ( 80A ) WRITE (6,610) BLANK(1:LOW), (DASH, I=1,MX), '+ ' RETURN ! Normal return C END! C C####################################################################### C REAL FUNCTION BICUBC (F, NBOX, DX, DY, DFDX, DFDY) C C Perform a type of bicubic interpolation in a grid of values. C For a point located DX, DY distant from the corner of the grid C (defined to be 1,1), return both the interpolated value and C its first derivatives with respect to x and y. C IMPLICIT NONE INTEGER NBOX REAL F(NBOX,NBOX), TEMP(4), DFDXT(4) C REAL DX, DY, DFDX, DFDY, C1, C2, C3, C4 INTEGER JY C C By construction, the point at which we want to estimate the function C will lie between the second and third columns, and between the second C and third rows of F, at a distance of (DX,DY) from the (2,2) element C of F. C DO JY=1,4 C1 = 0.5*(F(3,JY)-F(1,JY)) C4 = F(3,JY) - F(2,JY) - C1 C2 = 3.*C4 - 0.5*(F(4,JY)-F(2,JY)) + C1 C3 = C4 - C2 C4 = DX*C3 TEMP(JY) = DX*(DX*(C4+C2)+C1)+F(2,JY) DFDXT(JY)= DX*(C4*3.+2.*C2)+C1 END DO C1 = 0.5*(TEMP(3)-TEMP(1)) C4 = TEMP(3) - TEMP(2) - C1 C2 = 3.*C4 - 0.5*(TEMP(4)-TEMP(2)) + C1 C3 = C4 - C2 C4 = DY*C3 BICUBC = DY*(DY*(C4+C2)+C1)+TEMP(2) DFDY = DY*(C4*3.+2.*C2)+C1 C1 = 0.5*(DFDXT(3)-DFDXT(1)) C4 = DFDXT(3) - DFDXT(2) - C1 C2 = 3.*C4 - 0.5*(DFDXT(4)-DFDXT(2)) + C1 C3 = C4 - C2 DFDX = DY*(DY*(DY*C3+C2)+C1)+DFDXT(2) RETURN END! C C####################################################################### C REAL FUNCTION PROFIL (IPSTYP, DX, DY, PAR, DHDXC, DHDYC, . TERM, IDERIV) C C Compute the value of an ANALYTIC prfile for a point DX,DY distant C from the centroid. Return both the computed value and its C first derivatives with respect to x and y. If IDERIV .NE. 0, C return also the first derivatives with respect to all the parameters C defining the profile. C IMPLICIT NONE INTEGER MAXPAR, MAXPT PARAMETER (MAXPAR=6, MAXPT=4) C REAL PAR(MAXPAR), TERM(MAXPAR) REAL D(MAXPT,MAXPT), W(MAXPT,MAXPT) REAL X(MAXPT), XSQ(MAXPT), P1XSQ(MAXPT) C REAL EXP, DAOERF C REAL DX, DY, DHDXC, DHDYC, WFSQ, Y, WT, WF, ONEMP3 REAL RSQ, E, TALPHA, P1SQ, P2SQ, XY, DENOM REAL FUNC, YSQ, WP4FOD, P4FOD, F, P1P2, ERFX, DHDSX, ERFY REAL DEBY, DFBY, DBYX0, DBYY0 REAL DHDSY, ALPHA, P2YSQ INTEGER I, IPSTYP, IDERIV, IX, IY, NPT C DATA D / 0.00000000, 0.0, 0.0 , 0.0 , . -0.28867513, 0.28867513, 0.0 , 0.0 , . -0.38729833, 0.00000000, 0.38729833, 0.0 , . -0.43056816, -0.16999052, 0.16999052, 0.43056816/ DATA W / 1.00000000, 0.0 , 0.0 , 0.0 , . 0.50000000, 0.50000000, 0.0 , 0.0 , . 0.27777778, 0.44444444, 0.27777778, 0.0 , . 0.17392742, 0.32607258, 0.32607258, 0.17392742/ C PROFIL = 0. DHDXC = 0. DHDYC = 0. C IF (IDERIV .GT. 0) THEN DO I=1,MAXPAR TERM(I) = 0. END DO END IF C IF (IPSTYP .EQ. 1) THEN C C GAUSSIAN C C F = ERFX * ERFY / (PAR(1) * PAR(2)) C C PAR(1) is the HWHM in X; sigma(x) = 0.8493218 * HWHM C PAR(2) is the HWHM in Y; ditto C P1P2 = PAR(1)*PAR(2) ERFX = DAOERF(DX, 0., PAR(1), DHDXC, DHDSX) ERFY = DAOERF(DY, 0., PAR(2), DHDYC, DHDSY) PROFIL = ERFX*ERFY/P1P2 DHDXC = DHDXC*ERFY/P1P2 DHDYC = DHDYC*ERFX/P1P2 IF (IDERIV .GT. 0) THEN TERM(1) = (DHDSX-ERFX/PAR(1))*ERFY/P1P2 TERM(2) = (DHDSY-ERFY/PAR(2))*ERFX/P1P2 END IF ELSE IF (IPSTYP .EQ. 3) THEN C C MOFFAT FUNCTION BETA = 2.5 C BETA-1 C F = -------------------------------------------------------- C Ax * Ay * [1 + (X/Ax)**2 + (Y/Ay)**2 + (XY*Axy)]**BETA C C PAR(1) is the HWHM in x at y = 0: C C 1/2 = 1/[1 + (PAR(1)/Ax)**2]**BETA C so C 2**(1/BETA) - 1 = (PAR(1)/Ax)**2 C C Ax**2 = PAR(1)**2/[2**(1/BETA) - 1] C C When BETA = 2.5, Ax**2 = 3.129813 * PAR(1)**2 C C Hence, let us use C C 1 C F = --------------------------------------------------------------- C P(1)*P(2)*{1+0.3195079*[(X/P(1))**2+(Y/P(2))**2+(XY*P(3))]**2.5 C C neglecting a constant of proportionality. C ALPHA = 0.3195079 TALPHA = 0.6390158 ! 2.*ALPHA P1SQ = PAR(1)**2 P2SQ = PAR(2)**2 P1P2 = PAR(1)*PAR(2) XY = DX*DY C DENOM = 1. + ALPHA*(DX**2/P1SQ + DY**2/P2SQ + XY*PAR(3)) IF (DENOM .GT. 1.E4) RETURN FUNC = 1. / (P1P2 * DENOM**PAR(4)) IF (FUNC .GE. 0.046) THEN NPT = 4 ELSE IF (FUNC .GE. 0.0022) THEN NPT = 3 ELSE IF (FUNC .GE. 0.0001) THEN NPT = 2 ELSE IF (FUNC .GE. 1.E-10) THEN PROFIL = (PAR(4) - 1.) * FUNC P4FOD = PAR(4)*ALPHA*PROFIL/DENOM DHDXC = P4FOD*(2.*DX/P1SQ + DY*PAR(3)) DHDYC = P4FOD*(2.*DY/P2SQ + DX*PAR(3)) IF (IDERIV .GT. 0) THEN TERM(1) = (2.*P4FOD*DX**2/P1SQ-PROFIL)/PAR(1) TERM(2) = (2.*P4FOD*DY**2/P2SQ-PROFIL)/PAR(2) TERM(3) = - P4FOD*XY C TERM(4) = PROFIL*(1./(PAR(4)-1.)-ALOG(DENOM)) END IF RETURN ELSE RETURN END IF C DO IX=1,NPT X(IX) = DX+D(IX,NPT) XSQ(IX) = X(IX)**2 P1XSQ(IX) = XSQ(IX)/P1SQ END DO C DO IY=1,NPT Y = DY+D(IY,NPT) YSQ = Y**2 P2YSQ = YSQ/P2SQ DO IX=1,NPT WT = W(IY,NPT)*W(IX,NPT) XY=X(IX)*Y DENOM = 1. + ALPHA*(P1XSQ(IX) + P2YSQ + XY*PAR(3)) FUNC = (PAR(4) - 1.) / (P1P2 * DENOM**PAR(4)) P4FOD = PAR(4)*ALPHA*FUNC/DENOM WP4FOD = WT*P4FOD WF = WT*FUNC PROFIL = PROFIL + WF DHDXC = DHDXC + WP4FOD*(2.*X(IX)/P1SQ + Y*PAR(3)) DHDYC = DHDYC + WP4FOD*(2.*Y/P2SQ + X(IX)*PAR(3)) IF (IDERIV .GT. 0) THEN TERM(1) = TERM(1) + . (2.*WP4FOD*P1XSQ(IX)-WF)/PAR(1) TERM(2) = TERM(2) + . (2.*WP4FOD*P2YSQ-WF)/PAR(2) TERM(3) = TERM(3) - WP4FOD*XY C TERM(4) = TERM(4) + WF*(1./(PAR(4)-1.)-ALOG(DENOM)) END IF END DO END DO ELSE IF (IPSTYP .EQ. 5) THEN C C Penny function --- Gaussian core plus Lorentzian wings. The Lorentzian C is elongated along the x or y axis, the Gaussian may be tilted. C P1SQ = PAR(1)**2 P2SQ = PAR(2)**2 ONEMP3 = 1.-PAR(3) XY = DX*DY C RSQ = DX**2/P1SQ + DY**2/P2SQ IF (RSQ .GT. 1.E10) RETURN C F = 1./(1.+RSQ) RSQ = RSQ + XY*PAR(4) IF (RSQ .LT. 34.) THEN E = EXP(-0.6931472*RSQ) FUNC = PAR(3)*E + ONEMP3*F ELSE E = 0. FUNC = ONEMP3*F END IF C IF (FUNC .GE. 0.046) THEN NPT = 4 ELSE IF (FUNC .GE. 0.0022) THEN NPT = 3 ELSE IF (FUNC .GE. 0.0001) THEN NPT = 2 ELSE IF (FUNC .GE. 1.E-10) THEN PROFIL = FUNC DFBY = ONEMP3*F**2 DEBY = 0.6931472*PAR(3)*E DBYX0 = 2.*DX/P1SQ DBYY0 = 2.*DY/P2SQ DHDXC = DEBY*(DBYX0 + DY*PAR(4)) + DFBY*DBYX0 DHDYC = DEBY*(DBYY0 + DX*PAR(4)) + DFBY*DBYY0 IF (IDERIV .GT. 0) THEN DBYX0 = DBYX0*DX/PAR(1) DBYY0 = DBYY0*DY/PAR(2) DFBY = DFBY + DEBY TERM(1) = DFBY * DBYX0 TERM(2) = DFBY * DBYY0 TERM(3) = E - F TERM(4) = - DEBY * XY . / (0.5 - ABS(PAR(4))) END IF RETURN ELSE RETURN END IF C DO IX=1,NPT X(IX) = DX+D(IX,NPT) P1XSQ(IX) = X(IX)/P1SQ END DO C DO IY=1,NPT Y = DY+D(IY,NPT) P2YSQ = Y/P2SQ DO IX=1,NPT WT = W(IY,NPT)*W(IX,NPT) XY = X(IX)*Y RSQ = P1XSQ(IX)*X(IX) + P2YSQ*Y F = 1./(1.+RSQ) RSQ = RSQ + XY*PAR(4) IF (RSQ .LT. 34.) THEN E = EXP(-0.6931472*RSQ) FUNC = PAR(3)*E + ONEMP3*F DEBY = 0.6931472*WT*PAR(3)*E ELSE E = 0. FUNC = ONEMP3*F DEBY = 0. END IF PROFIL = PROFIL + WT*FUNC DFBY = WT*ONEMP3*F**2 DBYX0 = 2.*P1XSQ(IX) DBYY0 = 2.*P2YSQ DHDXC = DHDXC + . DEBY*(DBYX0 + DY*PAR(4)) + DFBY*DBYX0 DHDYC = DHDYC + . DEBY*(DBYY0 + DX*PAR(4)) + DFBY*DBYY0 IF (IDERIV .GT. 0) THEN DBYX0 = DBYX0*DX/PAR(1) DBYY0 = DBYY0*DY/PAR(2) TERM(1) = TERM(1) + (DFBY+DEBY)*DBYX0 TERM(2) = TERM(2) + (DFBY+DEBY)*DBYY0 TERM(3) = TERM(3) + WT*(E-F) TERM(4) = TERM(4) - DEBY * XY END IF END DO END DO ELSE IF (IPSTYP .EQ. 6) THEN C C Penny function --- Gaussian core plus Lorentzian wings. C The Lorentzian and Gaussian may be tilted in different C directions. C P1SQ = PAR(1)**2 P2SQ = PAR(2)**2 ONEMP3 = 1.-PAR(3) XY = DX*DY C RSQ = DX**2/P1SQ + DY**2/P2SQ DFBY = RSQ + PAR(5)*XY IF (DFBY .GT. 1.E10) RETURN F = 1./(1.+DFBY) C DEBY = RSQ + PAR(4)*XY IF (DEBY .LT. 34.) THEN E = EXP(-0.6931472*DEBY) ELSE E = 0. END IF C FUNC = PAR(3)*E + ONEMP3*F IF (FUNC .GE. 0.046) THEN NPT = 4 ELSE IF (FUNC .GE. 0.0022) THEN NPT = 3 ELSE IF (FUNC .GE. 0.0001) THEN NPT = 2 ELSE IF (FUNC .GE. 1.E-10) THEN PROFIL = FUNC DFBY = ONEMP3*F**2 DEBY = 0.6931472*PAR(3)*E DBYX0 = 2.*DX/P1SQ DBYY0 = 2.*DY/P2SQ DHDXC = DEBY*(DBYX0 + DY*PAR(4)) + . DFBY*(DBYX0 + DY*PAR(5)) DHDYC = DEBY*(DBYY0 + DX*PAR(4)) + . DFBY*(DBYY0 + DX*PAR(5)) IF (IDERIV .GT. 0) THEN DBYX0 = DBYX0*DX/PAR(1) DBYY0 = DBYY0*DY/PAR(2) TERM(5) = -DFBY * XY DFBY = DFBY + DEBY TERM(1) = DFBY * DBYX0 TERM(2) = DFBY * DBYY0 TERM(3) = E - F TERM(4) = - DEBY * XY END IF RETURN ELSE RETURN END IF C DO IX=1,NPT X(IX) = DX+D(IX,NPT) P1XSQ(IX) = X(IX)/P1SQ END DO C DO IY=1,NPT Y = DY+D(IY,NPT) P2YSQ = Y/P2SQ DO IX=1,NPT WT = W(IY,NPT)*W(IX,NPT) XY = X(IX)*Y RSQ = P1XSQ(IX)*X(IX) + P2YSQ*Y F = 1./(1.+RSQ+PAR(5)*XY) DEBY = RSQ + PAR(4)*XY IF (DEBY .LT. 34.) THEN E = EXP(-0.6931472*DEBY) FUNC = PAR(3)*E + ONEMP3*F DEBY = 0.6931472*WT*PAR(3)*E ELSE E = 0. FUNC = ONEMP3*F DEBY = 0. END IF PROFIL = PROFIL + WT*FUNC DFBY = WT*ONEMP3*F**2 DBYX0 = 2.*P1XSQ(IX) DBYY0 = 2.*P2YSQ DHDXC = DHDXC + . DEBY*(DBYX0 + DY*PAR(4)) + . DFBY*(DBYX0 + DY*PAR(5)) DHDYC = DHDYC + . DEBY*(DBYY0 + DX*PAR(4)) + . DFBY*(DBYY0 + DX*PAR(5)) IF (IDERIV .GT. 0) THEN DBYX0 = DBYX0*DX/PAR(1) DBYY0 = DBYY0*DY/PAR(2) TERM(1) = TERM(1) + (DFBY+DEBY)*DBYX0 TERM(2) = TERM(2) + (DFBY+DEBY)*DBYY0 TERM(3) = TERM(3) + WT*(E-F) TERM(4) = TERM(4) - DEBY * XY TERM(5) = TERM(5) - DFBY * XY END IF END DO END DO ELSE IF (IPSTYP .EQ. 2) THEN C C MOFFAT FUNCTION BETA = 1.5 C C BETA-1 C F = -------------------------------------------------------- C Ax * Ay * [1 + (X/Ax)**2 + (Y/Ay)**2 + (XY*Axy)]**BETA C C PAR(1) is the HWHM in x at y = 0: C C 1/2 = 1/[1 + (PAR(1)/Ax)**2]**BETA C so C 2**(1/BETA) - 1 = (PAR(1)/Ax)**2 C C Ax**2 = PAR(1)**2/[2**(1/BETA) - 1] C C When BETA = 1.5, Ax**2 = 1.7024144 * PAR(1)**2 C C Hence, let us use C C 1 C F = --------------------------------------------------------------- C P(1)*P(2)*{1+0.5874011*[(X/P(1))**2+(Y/P(2))**2+(XY*P(3))]**1.5 C C neglecting a constant of proportionality. C ALPHA = 0.5874011 TALPHA = 1.1748021 ! 2.*ALPHA P1SQ = PAR(1)**2 P2SQ = PAR(2)**2 P1P2 = PAR(1)*PAR(2) XY = DX*DY C DENOM = 1. + ALPHA*(DX**2/P1SQ + DY**2/P2SQ + XY*PAR(3)) IF (DENOM .GT. 5.E6) RETURN FUNC = 1. / (P1P2 * DENOM**PAR(4)) IF (FUNC .GE. 0.046) THEN NPT = 4 ELSE IF (FUNC .GE. 0.0022) THEN NPT = 3 ELSE IF (FUNC .GE. 0.0001) THEN NPT = 2 ELSE IF (FUNC .GE. 1.E-10) THEN PROFIL = (PAR(4) - 1.) * FUNC P4FOD = PAR(4)*ALPHA*PROFIL/DENOM DHDXC = P4FOD*(2.*DX/P1SQ + DY*PAR(3)) DHDYC = P4FOD*(2.*DY/P2SQ + DX*PAR(3)) IF (IDERIV .GT. 0) THEN TERM(1) = (2.*P4FOD*DX**2/P1SQ-PROFIL)/PAR(1) TERM(2) = (2.*P4FOD*DY**2/P2SQ-PROFIL)/PAR(2) TERM(3) = - P4FOD*XY C TERM(4) = PROFIL*(1./(PAR(4)-1.)-ALOG(DENOM)) END IF RETURN ELSE RETURN END IF C DO IX=1,NPT X(IX) = DX+D(IX,NPT) XSQ(IX) = X(IX)**2 P1XSQ(IX) = XSQ(IX)/P1SQ END DO C DO IY=1,NPT Y = DY+D(IY,NPT) YSQ = Y**2 P2YSQ = YSQ/P2SQ DO IX=1,NPT WT = W(IY,NPT)*W(IX,NPT) XY=X(IX)*Y DENOM = 1. + ALPHA*(P1XSQ(IX) + P2YSQ + XY*PAR(3)) FUNC = (PAR(4) - 1.)/ (P1P2 * DENOM**PAR(4)) P4FOD = PAR(4)*ALPHA*FUNC/DENOM WP4FOD = WT*P4FOD WF = WT*FUNC PROFIL = PROFIL + WF DHDXC = DHDXC + WP4FOD*(2.*X(IX)/P1SQ + Y*PAR(3)) DHDYC = DHDYC + WP4FOD*(2.*Y/P2SQ + X(IX)*PAR(3)) IF (IDERIV .GT. 0) THEN TERM(1) = TERM(1) + . (2.*WP4FOD*P1XSQ(IX)-WF)/PAR(1) TERM(2) = TERM(2) + . (2.*WP4FOD*P2YSQ-WF)/PAR(2) TERM(3) = TERM(3) - WP4FOD*XY C TERM(4) = TERM(4) + WF*(1./(PAR(4)-1.)-ALOG(DENOM)) END IF END DO END DO ELSE IF (IPSTYP .EQ. 4) THEN C C LORENTZ FUNCTION C 1 C F = -------------------------------------- C [1 + (X/Ax)**2 + (Y/Ay)**2 + (XY*Axy)] C C PAR(1) is the HWHM in x at y = 0. C P1SQ = PAR(1)**2 P2SQ = PAR(2)**2 P1P2 = PAR(1)*PAR(2) XY = DX*DY C DENOM = 1. + DX**2/P1SQ + DY**2/P2SQ + XY*PAR(3) IF (DENOM .GT. 1.E10) RETURN FUNC = 1. / DENOM IF (FUNC .GE. 0.046) THEN NPT = 4 ELSE IF (FUNC .GE. 0.0022) THEN NPT = 3 ELSE IF (FUNC .GE. 0.0001) THEN NPT = 2 ELSE IF (FUNC .GE. 1.E-10) THEN PROFIL = FUNC WFSQ = FUNC**2 DHDXC = WFSQ*(2.*DX/P1SQ + DY*PAR(3)) DHDYC = WFSQ*(2.*DY/P2SQ + DX*PAR(3)) IF (IDERIV .GT. 0) THEN TERM(1) = WFSQ*(2.*DX**2/P1SQ)/PAR(1) TERM(2) = WFSQ*(2.*DY**2/P2SQ)/PAR(2) TERM(3) = - WFSQ*XY END IF RETURN ELSE RETURN END IF C DO IX=1,NPT X(IX) = DX+D(IX,NPT) XSQ(IX) = X(IX)**2 P1XSQ(IX) = XSQ(IX)/P1SQ END DO C DO IY=1,NPT Y = DY+D(IY,NPT) YSQ = Y**2 P2YSQ = YSQ/P2SQ DO IX=1,NPT WT = W(IY,NPT)*W(IX,NPT) XY=X(IX)*Y DENOM = 1. + P1XSQ(IX) + P2YSQ + XY*PAR(3) FUNC = 1. / DENOM WF = WT*FUNC WFSQ = WF*FUNC PROFIL = PROFIL + WF DHDXC = DHDXC + WFSQ*(2.*X(IX)/P1SQ + Y*PAR(3)) DHDYC = DHDYC + WFSQ*(2.*Y/P2SQ + X(IX)*PAR(3)) IF (IDERIV .GT. 0) THEN TERM(1) = TERM(1) + WFSQ*(2.*P1XSQ(IX))/PAR(1) TERM(2) = TERM(2) + WFSQ*(2.*P2YSQ)/PAR(2) TERM(3) = TERM(3) - WFSQ*XY END IF END DO END DO ELSE write(*,*)'Invalid PSF type.' write(*,*)'Bye' stop END IF RETURN END! C C======================================================================= C FUNCTION NPARAM (IPSTYP, FWHM, LABEL, PAR, MAXPAR) CHARACTER*8 LABEL INTEGER MAXPAR REAL PAR(MAXPAR) PAR(1) = FWHM/2. PAR(2) = PAR(1) IF (IPSTYP .EQ. 1) THEN NPARAM = 2 LABEL = 'GAUSSIAN' ELSE IF (IPSTYP .EQ. 3) THEN NPARAM = 3 PAR(3) = 0. PAR(4) = 2.5 LABEL = 'MOFFAT25' ELSE IF (IPSTYP .EQ. 5) THEN NPARAM = 4 PAR(3) = 0.75 PAR(4) = 0.0 LABEL = 'PENNY1 ' ELSE IF (IPSTYP .EQ. 6) THEN NPARAM = 5 PAR(3) = 0.75 PAR(4) = 0.0 PAR(5) = 0.0 LABEL = 'PENNY2 ' ELSE IF (IPSTYP .EQ. 2) THEN NPARAM = 3 PAR(3) = 0. PAR(4) = 1.5 LABEL = 'MOFFAT15' ELSE IF (IPSTYP .EQ. 4) THEN NPARAM = 3 PAR(3) = 0. LABEL = 'LORENTZ ' ELSE write(*,*)'Invalid PSF type: '//CHAR(IPSTYP+48) END IF RETURN END! C C####################################################################### C REAL FUNCTION DAOERF (XIN, XO, BETA, DFDXO, DFDBET) C C Numerically integrate a Gaussian function C C F = EXP {-0.5*[(x-XO)/SIGMA]**2 }, C C from XIN-0.5 to XIN+0.5 using Gauss-Legendre integration. BETA C is the half-width at half-maximum, which is equal to 1.17741 * SIGMA. C Thus, C C F = EXP {-0.6931472*[(x-XO)/BETA]**2 }. C C Also: provide the first derivative of the integral with respect to C Xo and BETA. Use Gauss-Legendre integration. C C----------------------------------------------------------------------- C IMPLICIT NONE INTEGER MAXPT PARAMETER (MAXPT=4) C REAL DX(MAXPT,MAXPT), WT(MAXPT,MAXPT) C REAL EXP C REAL X, XSQ REAL XIN, XO, BETA, DFDXO, DFDBET, BETASQ, DELTAX, F, WF INTEGER NPT, I C DATA DX / 0.00000000, 0.0, 0.0 , 0.0 , . -0.28867513, 0.28867513, 0.0 , 0.0 , . -0.38729833, 0.00000000, 0.38729833, 0.0 , . -0.43056816, -0.16999052, 0.16999052, 0.43056816/ DATA WT / 1.00000000, 0.0 , 0.0 , 0.0 , . 0.50000000, 0.50000000, 0.0 , 0.0 , . 0.27777778, 0.44444444, 0.27777778, 0.0 , . 0.17392742, 0.32607258, 0.32607258, 0.17392742/ DAOERF = 0. DFDXO = 0. DFDBET = 0. BETASQ=BETA**2 DELTAX = XIN-XO C XSQ = DELTAX**2 F = XSQ/BETASQ IF (F .GT. 34.) RETURN F = EXP(-0.6931472*F) IF (F .GE. 0.046) THEN NPT = 4 ELSE IF (F .GE. 0.0022) THEN NPT = 3 ELSE IF (F .GE. 0.0001) THEN NPT = 2 ELSE IF (F .GE. 1.E-10) THEN DAOERF = F DFDXO = 1.3862944 * DELTAX * F / BETASQ DFDBET = 1.3862944 * XSQ * F / (BETASQ*BETA) RETURN ELSE RETURN END IF C DO I=1,NPT X = DELTAX + DX(I,NPT) XSQ = X**2 F = EXP(-0.6931472*XSQ/BETASQ) WF = WT(I,NPT)*F DAOERF = DAOERF+WF DFDXO = DFDXO + X*WF DFDBET = DFDBET + XSQ*WF END DO DFDXO = 1.3862944*DFDXO/BETASQ DFDBET = 1.3862944*DFDBET/(BETASQ*BETA) C RETURN END! C C####################################################################### C SUBROUTINE PKFIT (F, NX, NY, MAXBOX, X, Y, SCALE, SKY, RADIUS, . LOBAD, HIBAD, BRIGHT, IPSTYP, PAR, MAXPAR, NPAR, . PSF, MAXPSF, MAXEXP, NPSF, NEXP, NFRAC, . DELTAX, DELTAY, ERRMAG, CHI, SHARP, NITER, NCOL, NROW) C C======================================================================= C C This is the subroutine which does the actual one-star least-squares C profile fit for PEAK. C C OFFICIAL DAO VERSION: 1991 April 6 C C Arguments C F (INPUT) is an NX by NY array containing actual picture data. C C MAXBOX (INPUT) is the maximum value allowable for either NX or NY, C needed for the dimension statements below. PEAK and PSF will C provide different values of MAXBOX. C C SCALE (INPUT/OUTPUT) is the initial estimate of the brightness of C the star, expressed as a fraction of the brightness of the C PSF. Upon return, the final computed value of SCALE will C be passed back to the calling routine. C C X, Y (INPUT/OUTPUT) are the initial estimates of the centroid of C the star relative to the corner (1,1) of the subarray. Upon C return, the final computed values of X and Y will be passed C back to the calling routine. C C SKY (INPUT) is the local sky brightness value, carried on from C PHOTOMETRY via the data files. C C RADIUS (INPUT) is the fitting radius-- only pixels within RADIUS of C the instantaneous estimate of the star's centroid will be C included in the fit. C C LOBAD and HIBAD (INPUT) are bad pixel limits-- any pixel whose C brightness value falls outside this range will be presumed to C be bad, and will be ignored. C C BRIGHT (INPUT) contains the brightness normalization of the model C analytic point spread function C C PAR (INPUT) contains the values of the remaining NPARAM parameters C defining the analytic function which approximates the core of C the PSF. C C PSF (INPUT) is an NPSF by NPSF by NEXP+NFRAC look-up table C containing corrections from the analytic approximation of the C PSF to the true PSF. C C ERRMAG (OUTPUT) is the estimated standard error of the value of SCALE C returned by this routine. C C CHI (OUTPUT) is the estimated goodness-of-fit statistic: the ratio C of the observed pixel-to-pixel mean absolute deviation from C the profile fit, to the value expected on the basis of the C read-out noise and the photons/ADU (which are brought in C through COMMON block /ERROR/). C C SHARP (OUTPUT) is a goodness-of-fit statistic describing how much C broader the actual profile of the object appears than the C profile of the PSF. C C NITER (OUTPUT) is the number of iterations the solution required to C achieve convergence. If NITER = 50, the solution did not C converge. If for some reason a singular matrix occurs during C the least-squares solution, this will be flagged by setting C NITER = -1. C C======================================================================= C IMPLICIT NONE INTEGER MAXPSF, MAXEXP, MAXBOX, MAXPAR C C Parameter C C MAXPSF is the length of the side of the largest PSF look-up table C allowed. (Note: half-pixel grid spacing) C REAL C(3,3), V(3), CLAMP(3), DTOLD(3) REAL F(MAXBOX,MAXBOX), T(3), DT(3), NUMER REAL PSF(MAXPSF,MAXPSF,MAXEXP), PAR(MAXPAR) REAL USEPSF, AMAX1, AMIN1, ABS INTEGER MIN0, MAX0 C REAL LOBAD, DF, WT, DFDSIG, RHOSQ, RELERR, SIG, SIGSQ REAL FPOS, DVDXC, DVDYC, RSQ, DX, DXSQ, DY, DYSQ, DATUM REAL DENOM, SUMWT, CHIOLD, RADSQ, SHARP, CHI, ERRMAG REAL DELTAX, DELTAY, X, Y, PHPADU, RONOIS, PERR, PKERR REAL SCALE, SKY, RADIUS, HIBAD, BRIGHT INTEGER I, J, IX, IY, IXLO, IXHI, IYLO, IYHI, ISTAT INTEGER NPIX, NFRAC, NEXP, NITER, NPAR, NPSF, IPSTYP INTEGER NX, NY, NCOL, NROW LOGICAL CLIP, REDO C COMMON /ERROR/ PHPADU, RONOIS, PERR, PKERR C C----------------------------------------------------------------------- C C Initialize a few things for the solution. C RADSQ=RADIUS**2 DO 1010 I=1,3 CLAMP(I)=1. 1010 DTOLD(I)=0.0 CHIOLD=1. NITER=0 SHARP=0. CLIP = .FALSE. C C----------------------------------------------------------------------- C C Here begins the big least-squares loop. C 2000 NITER=NITER+1 C C Initialize things for this iteration. CHI and SHARP will be C goodness-of-fit indices. CHI will also be used in determining the C weights of the individual pixels. As the solution iterates, the new C value of CHI will be built up in the variable CHI, while a smoothed C value of CHI computed from the previous iteration will be carried C along in CHIOLD. C CHI=0.0 SUMWT=0.0 NUMER=0.0 DENOM=0.0 DO 2010 I=1,3 V(I)=0.0 ! Zero the vector of residuals DO 2010 J=1,3 2010 C(I,J)=0.0 ! Zero the normal matrix C C Choose the little box containing points inside the fitting radius. C IXLO=MAX0(1, INT(X-RADIUS)) IYLO=MAX0(1, INT(Y-RADIUS)) IXHI=MIN0(NX, INT(X+RADIUS)+1) IYHI=MIN0(NY, INT(Y+RADIUS)+1) C C Now build up the normal matrix and vector of residuals. C NPIX=0 DO 2090 IY=IYLO,IYHI DY=FLOAT(IY)-Y DYSQ=DY**2 DO 2090 IX=IXLO,IXHI DATUM=F(IX,IY) IF ((DATUM .LT. LOBAD) .OR. (DATUM .GT. HIBAD)) GO TO 2090 DX=FLOAT(IX)-X DXSQ=DX**2 C C DX and DY are the distance of this pixel from the centroid of the C star. Is this pixel truly inside the fitting radius? C RSQ=(DXSQ+DYSQ)/RADSQ IF (1.-RSQ .LE. 2.E-6) GO TO 2090 ! Prevents floating underflows C C The fitting equation is of the form C C Observed brightness = C [SCALE + delta(SCALE)] * [PSF + delta(Xcen)*d(PSF)/d(Xcen) + C delta(Ycen)*d(PSF)/d(Ycen) ] C C and is solved for the unknowns delta(SCALE) ( = the correction to C the brightness ratio between the program star and the PSF) and C delta(Xcen) and delta(Ycen) ( = corrections to the program star's C centroid). C C The point-spread function is equal to the sum of the integral under C a two-dimensional analytic profile plus values interpolated from C a look-up table. C T(1)=USEPSF(IPSTYP, DX, DY, BRIGHT, PAR, PSF, NPSF, NPAR, NEXP, . NFRAC, DELTAX, DELTAY, DVDXC, DVDYC) IF ((SCALE*T(1)+SKY .GT. HIBAD) .AND. (NITER.GE.3)) GO TO 2090 T(2)=SCALE*DVDXC T(3)=SCALE*DVDYC DF=F(IX,IY)-SKY-SCALE*T(1) C C DF is the residual of the brightness in this pixel from the PSF fit. C C The expected random error in the pixel is the quadratic sum of C the Poisson statistics, plus the readout noise, plus an estimated C error of PERERR% of the total brightness for the difficulty of flat- C fielding and bias-correcting the chip, plus an estimated error of C some fraction of the fourth derivative at the peak of the profile, C to account for the difficulty of accurately interpolating within the C point-spread function. The fourth derivative of the PSF is roughly C proportional to H/sigma**4 (sigma is the width parameter for C the stellar core); using the geometric mean of sigma(x) and sigma(y), C this becomes H/[sigma(x)*sigma(y)]**2. The ratio of the fitting C error to this quantity is estimated from a good-seeing CTIO frame to C be approximately 0.027. NOWADAYS, PAR(1) and PAR(2) are the HWHM, C not the sigma, so the coefficient is now of order 0.05 = 5.0%. C FPOS=AMAX1(0., F(IX,IY)-DF) C C FPOS = raw data minus residual = model-predicted value of the C intensity at this point (which presumably is non-negative). C SIGSQ=FPOS/PHPADU+RONOIS+(PERR*FPOS)**2+(PKERR*(FPOS-SKY))**2 SIG=SQRT(SIGSQ) RELERR=ABS(DF/SIG) C C SIG is the anticipated standard error of the intensity in this pixel, C including readout noise, Poisson photon statistics, and an estimate C of the standard error of interpolating within the PSF. C WT=5./(5.+RSQ/(1.-RSQ)) ! Weight as function of radius C C Now add this pixel into the quantities which go to make up the SHARP C index. C RHOSQ=DXSQ/PAR(1)**2+DYSQ/PAR(2)**2 C C Include in the sharpness index only those pixels within six C HWHMs of the centroid of the star. (This saves time and C floating underflows by excluding pixels which contribute less than C about one part in a million to the index.) C IF (RHOSQ .LE. 36.) THEN RHOSQ=0.6931472*RHOSQ DFDSIG=EXP(-RHOSQ)*(RHOSQ-1.) FPOS=AMAX1(0., F(IX,IY)-SKY)+SKY NUMER=NUMER+DFDSIG*DF/SIGSQ DENOM=DENOM+DFDSIG**2/SIGSQ END IF C C Derive the weight of this pixel. First of all, the weight depends C upon the distance of the pixel from the centroid of the star-- it C is determined from a function which is very nearly unity for radii C much smaller than the fitting radius, and which goes to zero for C radii very near the fitting radius. Then reject any pixels with C 10-sigma residuals (after the first iteration). C CHI=CHI+WT*RELERR SUMWT=SUMWT+WT C C Now the weight is scaled to the inverse square of the expected mean C error. C WT=WT/SIGSQ C C Reduce the weight of a bad pixel. A pixel having a residual of 2.5 C sigma gets reduced to half weight; a pixel having a residual of 5. C sigma gets weight 1/257. C IF (CLIP) WT=WT/(1.+(0.4*RELERR/CHIOLD)**8) C C Now add the pixel into the vector of residuals and the normal matrix. C DO 2030 I=1,3 V(I)=V(I)+DF*T(I)*WT DO 2030 J=1,3 2030 C(I,J)=C(I,J)+T(I)*T(J)*WT C NPIX=NPIX+1 2090 CONTINUE ! End of loop over pixels C C Compute the (robust) goodness-of-fit index CHI. C IF (SUMWT .GT. 3.) THEN CHI=1.2533141*CHI/SQRT(SUMWT*(SUMWT-3.)) C C CHI is pulled toward its expected value of unity before being stored C in CHIOLD to keep the statistics of a small number of pixels from C completely dominating the error analysis. C CHIOLD=((SUMWT-3.)*CHI+3.)/SUMWT ELSE CHI = 1. CHIOLD = 1. END IF C C Compute the parameter corrections and check for convergence. C IF (NPIX .LT. 3) THEN NITER=-1 RETURN END IF CALL INVERS (C, 3, 3, ISTAT) IF (ISTAT .EQ. 0) GO TO 2100 NITER=-1 RETURN C C Everything OK so far. C 2100 CALL VMUL (C, 3, 3, V, DT) C C In the beginning, the brightness of the star will not be permitted C to change by more than two magnitudes per iteration (that is to say, C if the estimate is getting brighter, it may not get brighter by C more than 525% per iteration, and if it is getting fainter, it may C not get fainter by more than 84% per iteration). The x and y C coordinates of the centroid will be allowed to change by no more C than one-half pixel per iteration. Any time that a parameter C correction changes sign, the maximum permissible change in that C parameter will be reduced by a factor of 2. C DO 2110 I=1,3 IF (DTOLD(I)*DT(I) .LT. 0.) THEN CLAMP(I)=0.5*CLAMP(I) ELSE CLAMP(I)=MIN(1., 1.1*CLAMP(I)) END IF 2110 DTOLD(I)=DT(I) C SCALE=SCALE+DT(1)/ . (1.+AMAX1(DT(1)/(5.25*SCALE),-DT(1)/(0.84*SCALE))/CLAMP(1)) X=X+DT(2)/(1.+ABS(DT(2))/(0.4*CLAMP(2))) Y=Y+DT(3)/(1.+ABS(DT(3))/(0.4*CLAMP(3))) IF (NITER .LE. 1) GO TO 2000 REDO=.FALSE. C C Convergence criteria: if the most recent computed correction to the C brightness is larger than 0.01% or than 0.05 * sigma(brightness), C whichever is larger, OR if the absolute change in X or Y is C greater than 0.001 pixels, convergence has not been achieved. C ERRMAG=CHIOLD*SQRT(C(1,1)) IF (CLIP) THEN IF (ABS(DT(1)) .GT. . AMAX1( 0.05*ERRMAG, 0.0001*SCALE )) THEN REDO=.TRUE. ELSE IF (AMAX1(ABS(DT(2)),ABS(DT(3))) .GT. 0.001) THEN REDO=.TRUE. END IF ELSE IF (ABS(DT(1)) .GT. . AMAX1( ERRMAG, 0.002*SCALE )) THEN REDO = .TRUE. ELSE IF (AMAX1(ABS(DT(2)),ABS(DT(3))) .GT. 0.02) THEN REDO = .TRUE. END IF END IF C IF (NITER .GE. 50) REDO = .FALSE. IF (REDO) GO TO 2000 IF ((NITER .LT. 50) .AND. (.NOT. CLIP)) THEN CLIP = .TRUE. DO I=1,3 DTOLD(I) = 0. CLAMP(I) = AMAX1(CLAMP(I), 0.25) END DO GO TO 2000 ELSE SHARP=1.4427*PAR(1)*PAR(2)*NUMER/(BRIGHT*SCALE*DENOM) SHARP=AMIN1(99.999,AMAX1(SHARP,-99.999)) RETURN END IF C END! c subroutine ovrwrt (line, iwhich) character*(*) line if (iwhich .eq. 1) then write (6,1) line 1 format (/a) else if (iwhich .eq. 2) then write (6,2) line, char(13) 2 format (a, a1, $) else if (iwhich .eq. 3) then write (6,3) line 3 format (a) else write (6,4) line, char(13) 4 format (/a, a1, $) end if call flush (6) return end C C####################################################################### C REAL FUNCTION USEPSF (IPSTYP, DX, DY, BRIGHT, PAR, PSF, . NPSF, NPAR, NEXP, NFRAC, DELTAX, DELTAY, DVDXC, DVDYC) C C Evaluate the PSF for a point distant DX, DY from the center of a C star located at relative frame coordinates DELTAX, DELTAY. C IMPLICIT NONE INTEGER MAXPSF, MAXPAR, MAXEXP PARAMETER (MAXPSF=207, MAXPAR=6, MAXEXP=10) C REAL PAR(*), PSF(MAXPSF,MAXPSF,*), JUNK(MAXEXP) C REAL PROFIL, BICUBC C REAL MIDDLE, BRIGHT, DVDXC, DVDYC, DELTAX, DELTAY, XX, YY REAL DX, DY, CORR, DFDX, DFDY INTEGER K, LX, LY, IPSTYP INTEGER NFRAC, NTERM, NPSF, NEXP, NPAR C NTERM = NEXP + NFRAC USEPSF = BRIGHT*PROFIL(IPSTYP, DX, DY, PAR, DVDXC, DVDYC, . JUNK, 0) DVDXC = BRIGHT*DVDXC DVDYC = BRIGHT*DVDYC IF (NTERM .LT. 0) RETURN MIDDLE = (NPSF+1)/2 C C The PSF look-up tables are centered at (MIDDLE, MIDDLE). C IF (NEXP .GE. 0) THEN JUNK(1) = 1. IF (NEXP .GE. 2) THEN JUNK(2) = DELTAX JUNK(3) = DELTAY IF (NEXP .GE. 4) THEN JUNK(4) = 1.5*DELTAX**2-0.5 JUNK(5) = DELTAX*DELTAY JUNK(6) = 1.5*DELTAY**2-0.5 IF (NEXP .GE. 7) THEN JUNK(7) = DELTAX*(5.*JUNK(4)-2.)/3. JUNK(8) = JUNK(4)*DELTAY JUNK(9) = DELTAX*JUNK(6) JUNK(10) = DELTAY*(5.*JUNK(6)-2.)/3. END IF END IF END IF END IF C C IF (NFRAC .GT. 0) THEN C J = NEXP+1 C JUNK(J) = -2.*(DX - REAL(NINT(DX))) C J = J+1 C JUNK(J) = -2.*(DY - REAL(NINT(DY))) C J = J+1 C JUNK(J) = 1.5*JUNK(J-2)**2 - 0.5 C J = J+1 C JUNK(J) = JUNK(J-3)*JUNK(J-2) C J = J+1 C JUNK(J) = 1.5*JUNK(J-3)**2 - 0.5 C END IF XX = 2.*DX+MIDDLE LX = INT(XX) YY = 2.*DY+MIDDLE LY = INT(YY) C C This point in the stellar profile lies between columns LX and LX+1, C and between rows LY and LY+1 in the look-up tables. C DO K=1,NTERM CORR = BICUBC(PSF(LX-1,LY-1,K), MAXPSF, . XX-REAL(LX), YY-REAL(LY), DFDX, DFDY) USEPSF = USEPSF + JUNK(K)*CORR DVDXC = DVDXC-JUNK(K)*DFDX DVDYC = DVDYC-JUNK(K)*DFDY END DO RETURN END! C function indexpos(head,f1) character*80 head(72),f1*8 do 10 indexpos=1,72 10 if(head(indexpos)(1:8).eq.f1)return end