parameter (n2048=4064) PARAMETER (MAXEXP=10, MAXPSF=207, MAXPAR=6) real par(MAXPAR) real a_b(n2048*n2048) character*60 f1,f2,f3 character*80 head(36) INTEGER NCOL, NROW logical logi if(iargc().lt.1)then write(*,*) write(*,*) write(*,*)' ********** DAO_SUBSTAR 2000.7 ***********' write(*,*) write(*,*) write(*,*)' Usage: dao_substar file.fit' write(*,*) write(*,*) stop endif call getarg(1,f1) 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 !' k=index(f1,'.') f2=f1(1:k)//'psf' ! 1 read *.psf open(1,file=f2,status='old') f2=f1(1:k)//'ap' ! 2 read *.ap open(2,file=f2,status='old') f2=f1(1:k)//'lst' ! 3 read *.lst open(3,file=f2,status='old') call readfits(f1,a_b,ncol,nrow) C wtach=-2. call SUBSTR(PAR,MAXPAR,MAXPSF,MAXEXP,a_b,NCOL,NROW,WATCH) C f3=f1(1:lnblnk(f1))//'s'//char(0) open(1,file=f1,status='UNKNOWN',access='direct',recl=2880) i=0 502 i=i+1 read(1,rec=i)head call writefits(f3,head,2880,i) if(indexpos(head,"END ",36).gt.36)goto 502 close(1) call swap4(a_b,ncol*nrow*4) call writefits(f3,a_b,ncol*nrow*4,0) write(*,*)' ok!' END! C SUBROUTINE SUBSTR (PAR, MAXPAR, MAXPSF, MAXEXP, . F, NCOL, NROW, WATCH) C C======================================================================= C C This subroutine scales and shifts the point-spread function according C to each star's magnitude and centroid, and subtracts the resulting C profile from a copy of the original picture. C C OFFICIAL DAO VERSION: 1991 April 18 C C Arguments C C WATCH (INPUT) governs whether information relating to the progress C of the reductions is to be typed on the terminal screen C during execution. C C======================================================================= C IMPLICIT NONE INTEGER MAXBOX, MAXPSF, MAXPAR, MAXEXC, MAXEXP PARAMETER (MAXBOX=69, MAXEXC=300) C C Parameters C C MAXBOX is the side of the square subarray containing the largest C (circular) PSF that can be subtracted from the picture. C C MAXPSF is the largest permissible number of elements on a side of the C (square) look-up table for the point-spread function. C C MAXBOX = (MAXPSF-7)/2. C INTEGER NCOL, NROW REAL PAR(MAXPAR) REAL F(NCOL,NROW), PSF(MAXPSF,MAXPSF,MAXEXP) INTEGER IEXC(MAXEXC) C REAL USEPSF INTEGER RDPSF, MAX0, MIN0, INT C CHARACTER*80 LINE CHARACTER*30 SUBPIC CHARACTER*30 COOFIL, MAGFIL, PSFFIL, PROFIL, GRPFIL, SWITCH CHARACTER CASE*5, ANSWER*1 REAL LOBAD, X, Y, STRMAG, SKY, DX, DY, DYSQ, SCALE, COL, ROW REAL DELTAX, DELTAY, DVDXC, DVDYC, XPSF, YPSF, PSFRAD, DUM REAL HIBAD, THRESH, AP1, PHPADU, READNS, FRAD, WATCH, PSFMAG REAL BRIGHT, PSFRSQ INTEGER I, J, L, NL, IDUM, NEXC, LX, LY, NX, NY, ISTAR, ID INTEGER ISTAT, IPSTYP, NPSF, NPAR, NEXP, NFRAC C COMMON /FILNAM/ COOFIL, MAGFIL, PSFFIL, PROFIL, GRPFIL C C----------------------------------------------------------------------- C C SECTION 1 C C Get file names and set up the needed numerical constants. C c CALL TBLANK ! Type a blank line c 950 CALL GETNAM ('File with the PSF:', PSFFIL) c IF ((PSFFIL .EQ. 'END OF FILE') .OR. c . (PSFFIL .EQ. 'GIVE UP')) THEN c PSFFIL = ' ' c RETURN c END IF ISTAT = RDPSF (1, IPSTYP, PAR, MAXPAR, NPAR, . PSF, MAXPSF, MAXEXP, NPSF, NEXP, NFRAC, . PSFMAG, BRIGHT, XPSF, YPSF) c IF (ISTAT .LT. 0) THEN c PSFFIL = 'GIVE UP' c GO TO 950 c END IF C PSFRAD=(REAL(NPSF-1)/2. - 1.)/2. PSFRSQ=PSFRAD**2 COL=FLOAT(NCOL) ROW=FLOAT(NROW) C c CALL CLFILE (2) close(11) C c 1015 CALL GETNAM ('File with photometry:', PROFIL) c IF ((PROFIL .EQ. 'END OF FILE') .OR. c . (PROFIL .EQ. 'GIVE UP')) THEN c PROFIL = ' ' c RETURN c END IF c CALL INFILE (2, PROFIL, ISTAT) c IF (ISTAT .NE. 0) THEN c CALL STUPID ('Error opening input file '//PROFIL) c PROFIL = 'GIVE UP' c GO TO 1015 c END IF C c CALL GETYN ('Do you have stars to leave in?', ANSWER) c IF (ANSWER .EQ. 'Y') THEN c SUBPIC = SWITCH(PROFIL, CASE('.lst')) c 1025 CALL GETNAM ('File with star list:', SUBPIC) c IF ((SUBPIC .EQ. 'END OF FILE') .OR. (SUBPIC .EQ. c . 'GIVE UP')) THEN c CALL CLFILE (2) c RETURN c END IF c CALL INFILE (3, SUBPIC, ISTAT) c IF (ISTAT .NE. 0) THEN c CALL STUPID ('Error opening file '//SUBPIC) c SUBPIC = 'GIVE UP' c GO TO 1025 c END IF c c CALL CHECK (3, NL) ANSWER = 'Y' NEXC = 0 1035 NEXC = NEXC+1 if(NEXC.GT.MAXEXC)goto 1036 read(3,*,end=1036)IEXC(NEXC) goto 1035 1036 close(3) NEXC=NEXC-1 write(*,*)NEXC c NEXC = 0 c 1035 NEXC = NEXC+1 c 1036 CALL RDSTAR (3, NL, IEXC(NEXC), DUM, DUM, DUM, DUM) c IF (IEXC(NEXC) .EQ. 0) GO TO 1036 c IF (IEXC(NEXC) .LT. 0) THEN c NEXC = NEXC-1 c ELSE c IF (NEXC .LT. MAXEXC) GO TO 1035 c END IF c close(3) c write(*,*)NEXC,MAXEXC c IF (NEXC .LE. 0) ANSWER = 'N' c CALL CLFILE (3) c END IF C c CALL RDHEAD (2, NL, IDUM, IDUM, LOBAD, HIBAD, THRESH, AP1, c . PHPADU, READNS, FRAD) read(2,'(a)')i read(2,*)NL,IDUM,IDUM,LOBAD,HIBAD,THRESH,AP1,PHPADU,READNS,FRAD C c SUBPIC=SWITCH(PROFIL, CASE('s')) c CALL GETNAM ('Name for subtracted image:', SUBPIC) c IF (SUBPIC .EQ. 'END OF FILE') GO TO 9010 ! CTRL-Z was entered C C Copy the input picture verbatim into the output picture. C C CALL COPPIC (SUBPIC, F, NCOL, NROW, ISTAT) c CALL COPPIC (SUBPIC, ISTAT) c IF (ISTAT .NE. 0) THEN c CALL STUPID ('Error creating output image.') c CALL CLPIC ('COPY') c CALL CLFILE (2) c RETURN c END IF C IF (WATCH .GT. 0.5) THEN c CALL TBLANK write(*,*) CALL OVRWRT (' Star', 1) END IF C C Read the entire image in from the copy. C c LX = 1 c LY = 1 c NX = NCOL c NY = NROW c CALL RDARAY ('COPY', LX, LY, NX, NY, NCOL, F, ISTAT) C C----------------------------------------------------------------------- C C SECTION 2 C C Loop over stars. C ISTAR=0 2000 ISTAR=ISTAR+1 2010 CALL RDSTAR (2, NL, ID, X, Y, STRMAG, SKY) IF (ID .LT. 0) GO TO 9000 ! End-of-file encountered IF (ID .EQ. 0) GO TO 2010 ! Ignore a blank line IF ((1.-X .GT. PSFRAD) .OR. (X-COL .GT. PSFRAD) .OR. . (1.-Y .GT. PSFRAD) .OR. (Y-ROW .GT. PSFRAD)) GO TO 2000 IF (ANSWER .EQ. 'Y') THEN L = NEXC DO I=1,L IF (ID .EQ. IEXC(I)) THEN IEXC(I) = IEXC(NEXC) NEXC = NEXC-1 IF (NEXC .EQ. 0) ANSWER = 'N' GO TO 2010 END IF END DO END IF IF (STRMAG .GE. 99.) GO TO 2000 ! Ignore a bad star IF ((WATCH .GT. 0.5) .AND. (MOD(ISTAR,20) .EQ. 0)) THEN WRITE (LINE,620) ISTAR 620 FORMAT (I5) CALL OVRWRT (LINE(1:5), 2) END IF DELTAX=(X-1.)/XPSF-1. DELTAY=(Y-1.)/YPSF-1. LX=MAX0(1, INT(X-PSFRAD)+1) LY=MAX0(1, INT(Y-PSFRAD)+1) NX=MIN0(NCOL, INT(X+PSFRAD)) NY=MIN0(NROW, INT(Y+PSFRAD)) SCALE=10.**(0.4*(PSFMAG-STRMAG)) C C Subtract the shifted scaled PSF C DO 2030 J=LY,NY DY=FLOAT(J)-Y DYSQ=DY**2 DO 2020 I=LX,NX IF (F(I,J) .GT. HIBAD) GO TO 2020 DX=FLOAT(I)-X IF (DX**2+DYSQ .GE. PSFRSQ) THEN IF (DX .GT. 0.) GO TO 2030 ELSE F(I,J)=F(I,J)-SCALE*USEPSF(IPSTYP, DX, DY, BRIGHT, PAR, . PSF, NPSF, NPAR, NEXP, NFRAC, DELTAX, DELTAY, . DVDXC, DVDYC) END IF 2020 CONTINUE 2030 CONTINUE GO TO 2000 ! Go to next star C C----------------------------------------------------------------------- C C Normal return. C 9000 CONTINUE IF (WATCH .GT. 0.5) THEN WRITE (LINE,620) ISTAR-1 CALL OVRWRT (LINE(1:5), 2) END IF IF (WATCH .GT. 0.5) CALL OVRWRT (' ', 4) C C Write the modified image back into the copy. C c LX = 1 c LY = 1 c NX = NCOL c NY = NROW c CALL WRARAY ('COPY', LX, LY, NX, NY, NCOL, F, ISTAT) c CALL CLPIC ('COPY') c IF (WATCH .GT. -1.5) write(*,*)' done. ' 9010 close(2) RETURN C 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 QUICK (DATUM, N, INDEX) C C======================================================================= C C A quick-sorting algorithm suggested by the discussion on pages 114-119 C of THE ART OF COMPUTER PROGRAMMING, Vol. 3, SORTING AND SEARCHING, by C D.E. Knuth, which was referenced in Don Wells' subroutine QUIK. This C is my own attempt at encoding a quicksort-- PBS. C C Arguments C C DATUM (INPUT/OUTPUT) is a vector of dimension N containing randomly C ordered real data upon input. Upon output the elements of C DATUM will be in order of increasing value. C C C INDEX (OUTPUT) is an integer vector of dimension N. Upon return to C the calling program the i-th element of INDEX will tell where C the i-th element of the sorted vector DATUM had been BEFORE C DATUM was sorted. C C======================================================================= C IMPLICIT NONE INTEGER MAXSTK, N PARAMETER (MAXSTK=28) C C Parameter C C MAXSTK is the maximum number of entries the stack can contain. C A limiting stack length of 14 restricts this quicksort C subroutine to vectors of maximum length of order 32,768 C (= 2**15). REAL DATUM(N) INTEGER INDEX(N), STKLO(MAXSTK), STKHI(MAXSTK) C REAL DKEY INTEGER I, HI, LO, NSTAK, LIMLO, LIMHI, IKEY C C Initialize INDEX. C DO I=1,N INDEX(I)=I END DO C C Initialize the pointers. C NSTAK=0 LIMLO=1 LIMHI=N C 100 DKEY=DATUM(LIMLO) IKEY=INDEX(LIMLO) C TYPE *, 'LO =', LIMLO, ' HI =', LIMHI C C Compare all elements in the sub-vector between LIMLO and LIMHI with C the current key datum. C LO=LIMLO HI=LIMHI 101 CONTINUE C IF (LO .EQ. HI)GO TO 200 C IF (DATUM(HI) .LE. DKEY) GO TO 109 HI=HI-1 C C The pointer HI is to be left pointing at a datum SMALLER than the C key, which is intended to be overwritten. C GO TO 101 C 109 DATUM(LO)=DATUM(HI) INDEX(LO)=INDEX(HI) LO=LO+1 110 CONTINUE C IF (LO .EQ. HI) GO TO 200 C IF (DATUM(LO) .GE. DKEY) GO TO 119 C LO=LO+1 GO TO 110 C 119 DATUM(HI)=DATUM(LO) INDEX(HI)=INDEX(LO) HI=HI-1 C C The pointer LO is to be left pointing at a datum LARGER than the C key, which is intended to be overwritten. C GO TO 101 C 200 CONTINUE C C LO and HI are equal, and point at a value which is intended to C be overwritten. Since all values below this point are less than C the key and all values above this point are greater than the key, C this is where we stick the key back into the vector. C DATUM(LO)=DKEY INDEX(LO)=IKEY C DO 1666 I=LIMLO,LO-1 C1666 TYPE *, DATUM(I) C TYPE *, DATUM(LO), ' KEY' C DO 2666 I=LO+1,LIMHI C2666 TYPE *, DATUM(I) C C At this point in the subroutine, all data between LIMLO and LO-1, C inclusive, are less than DATUM(LO), and all data between LO+1 and C LIMHI are larger than DATUM(LO). C C If both subarrays contain no more than one element, then take the most C recent interval from the stack (if the stack is empty, we're done). C If the larger of the two subarrays contains more than one element, and C if the shorter subarray contains one or no elements, then forget the C shorter one and reduce the other subarray. If the shorter subarray C contains two or more elements, then place the larger subarray on the C stack and process the subarray. C IF (LIMHI-LO .GT. LO-LIMLO) GO TO 300 C C Case 1: the lower subarray is longer. If it contains one or no C elements then take the most recent interval from the stack and go C back and operate on it. C IF (LO-LIMLO .LE. 1) GO TO 400 C C If the upper (shorter) subinterval contains one or no elements, then C process the lower (longer) one, but if the upper subinterval contains C more than one element, then place the lower (longer) subinterval on C the stack and process the upper one. C IF (LIMHI-LO .GE. 2) GO TO 250 C C Case 1a: the upper (shorter) subinterval contains no or one elements, C so we go back and operate on the lower (longer) subinterval. C LIMHI=LO-1 GO TO 100 C 250 CONTINUE C C Case 1b: the upper (shorter) subinterval contains at least two C elements, so we place the lower (longer) subinterval on the stack and C then go back and operate on the upper subinterval. C NSTAK=NSTAK+1 IF (NSTAK .GT. MAXSTK) THEN write(*,*)'Stack overflow in QUICK. Increase MAXSTK.' END IF STKLO(NSTAK)=LIMLO STKHI(NSTAK)=LO-1 LIMLO=LO+1 GO TO 100 300 CONTINUE C C Case 2: the upper subarray is longer. If it contains one or no C elements then take the most recent interval from the stack and C operate on it. C IF (LIMHI-LO .LE. 1) GO TO 400 C C If the lower (shorter) subinterval contains one or no elements, then C process the upper (longer) one, but if the lower subinterval contains C more than one element, then place the upper (longer) subinterval on C the stack and process the lower one. C IF (LO-LIMLO .GE. 2) GO TO 350 C C Case 2a: the lower (shorter) subinterval contains no or one elements, C so we go back and operate on the upper (longer) subinterval. C LIMLO=LO+1 GO TO 100 C 350 CONTINUE C C Case 2b: the lower (shorter) subinterval contains at least two C elements, so we place the upper (longer) subinterval on the stack and C then go back and operate on the lower subinterval. C NSTAK=NSTAK+1 IF (NSTAK .GT. MAXSTK) THEN write(*,*)'Stack overflow in QUICK. Increase MAXSTK.' END IF STKLO(NSTAK)=LO+1 STKHI(NSTAK)=LIMHI LIMHI=LO-1 GO TO 100 C 400 CONTINUE C C Take the most recent interval from the stack. If the stack happens C to be empty, we are done. C IF (NSTAK .LE. 0) THEN RETURN ! Normal return END IF C TYPE *, 'POP: ', NSTAK, STKLO(NSTAK), STKHI(NSTAK) LIMLO=STKLO(NSTAK) LIMHI=STKHI(NSTAK) NSTAK=NSTAK-1 GO TO 100 C END! C C####################################################################### C SUBROUTINE RECTFY (X, NSTAR, INDEX, HOLD) IMPLICIT NONE C REAL X(*), HOLD(*) INTEGER INDEX(*) C INTEGER I, NSTAR C DO I=1,NSTAR HOLD(I)=X(I) END DO DO I=1,NSTAR X(I)=HOLD(INDEX(I)) END DO RETURN END! C SUBROUTINE RECTFYI (ID, NSTAR, INDEX, HOLD) IMPLICIT NONE C REAL HOLD(*) INTEGER INDEX(*),ID(*) C INTEGER I, NSTAR C DO I=1,NSTAR HOLD(I)=ID(I) END DO DO I=1,NSTAR ID(I)=HOLD(INDEX(I)) END DO RETURN END! C C####################################################################### C REAL FUNCTION PCTILE(DATUM,N,NPCT) C C======================================================================= C C This is a modification of a quick-sorting algorithm, which is intended C to take in a vector of numbers, and return the value of the PCT-th C percentile in that vector: C C DATUM (input real vector) containing real data. C N (input integer) number of elements in DATUM. C NPCT (input integer) element of sorted vector whose value C is desired. C PCTILE (output real) the value of the NPCT-th element C in the sorted vector DATUM. C C----------------------------------------------------------------------- C C The quick-sorting algorithm was suggested by the discussion on pages C 114-119 of THE ART OF COMPUTER PROGRAMMING, Vol. 3, SORTING AND C SEARCHING, by D.E. Knuth, which was referenced in Don Wells' C subroutine QUIK. This is my own attempt at encoding a quicksort-- C PBS. C C The array DATUM contains randomly ordered data. C IMPLICIT NONE REAL DATUM(*) C INTEGER MIN0, MAX0 C REAL DKEY INTEGER LO, HI, N, NPCT, LIMLO, LIMHI, NTARG C C Which element of the sorted array will we be interested in? C NTARG=MAX0(1,MIN0(N,NPCT)) C C Initialize the pointers. C LIMLO=1 LIMHI=N C 100 DKEY=DATUM(LIMLO) C TYPE *,'LOW=',LIMLO,' HIGH=',LIMHI,' KEY=',DKEY C C Compare all elements in the sub-vector between LIMLO and LIMHI with C the current key datum. C LO=LIMLO HI=LIMHI 101 CONTINUE C C If LO equals HI, we have tested all the elements in the current search C interval. C IF(LO.EQ.HI)GO TO 200 IF(DATUM(HI).LE.DKEY)GO TO 109 HI=HI-1 C C The pointer HI is to be left pointing at a datum SMALLER than the C key, which is intended to be overwritten. C GO TO 101 C 109 DATUM(LO)=DATUM(HI) LO=LO+1 110 CONTINUE IF(LO.EQ.HI)GO TO 200 IF(DATUM(LO).GE.DKEY)GO TO 119 LO=LO+1 GO TO 110 C 119 DATUM(HI)=DATUM(LO) HI=HI-1 C C The pointer LO is to be left pointing at a datum LARGER than the C key, which is intended to be overwritten. C GO TO 101 C 200 CONTINUE C C LO and HI are equal, and point at a value which is intended to C be overwritten. Since all values below this point are less than C the key and all values above this point are greater than the key, C this is where we stick the key back into the vector. C DATUM(LO)=DKEY C DO 1666 I=LIMLO,LO-1 C1666 TYPE *,DATUM(I) C TYPE *,DATUM(LO),' KEY' C DO 2666 I=LO+1,LIMHI C2666 TYPE *,DATUM(I) C C At this point in the subroutine, all data between LIMLO and LO-1, C inclusive, are less than DATUM(LO), and all data between LO+1 and C LIMHI are larger than DATUM(LO). If LO = NTARG, then DATUM(LO) is C the value we are looking for. If NTARG < LO, then we want to sort the C values of DATUM from LIMLO to LO-1, inclusive, whereas if NTARG > LO, C then we want to sort the values of DATUM from LO+1 to LIMHI, C inclusive. C C TYPE *,'NTARG=',NTARG,' LO=',LO IF(NTARG-LO)300,900,400 300 LIMHI=LO-1 GO TO 100 400 LIMLO=LO+1 GO TO 100 900 PCTILE=DATUM(LO) RETURN END! C C####################################################################### C REAL FUNCTION SMLLST (DATUM, N, M) IMPLICIT NONE REAL DATUM(*) INTEGER I, J, N, M, LEAST DO J=1,M SMLLST = DATUM(J) LEAST = J DO I=J,N IF (DATUM(I) .LT. SMLLST) THEN SMLLST = DATUM(I) LEAST = I END IF END DO DATUM(LEAST) = DATUM(J) DATUM(J) = SMLLST END DO RETURN END! C C####################################################################### C REAL FUNCTION BIGGST (DATUM, N, M) IMPLICIT NONE REAL DATUM(*) INTEGER I, J, N, M, MOST DO J=N,N-M+1,-1 MOST = J BIGGST = DATUM(MOST) DO I=1,J IF (DATUM(I) .GT. BIGGST) THEN BIGGST = DATUM(I) MOST = I END IF END DO DATUM(MOST) = DATUM(J) DATUM(J) = BIGGST END DO RETURN END! C C####################################################################### C INTEGER FUNCTION RDPSF (id, IPSTYP, PAR, MAXPAR, NPAR, . PSF, MAXPSF, MAXEXP, NPSF, NEXP, NFRAC, . PSFMAG, BRIGHT, XPSF, YPSF) C C Read in the point-spread function C IMPLICIT NONE integer id INTEGER MAXPSF, MAXTYP, MAXPAR, MAXEXP PARAMETER (MAXTYP=6) C REAL PAR(MAXPAR), PSF(MAXPSF,MAXPSF,MAXEXP) C CHARACTER*8 LABEL, CHECK REAL PSFMAG, BRIGHT, XPSF, YPSF INTEGER I, J, K, IPSTYP, NPSF, NPAR, NEXP, NFRAC, ISTAT INTEGER NTERM, NPARAM C READ(id,302,IOSTAT=ISTAT)LABEL,NPSF,NPAR,NEXP,NFRAC,PSFMAG, . BRIGHT, XPSF, YPSF 302 FORMAT (1X, A8, 4I5, F9.3, F15.3, 2F9.1) IF (ISTAT .NE. 0) THEN write(*,*)'Error reading PSF.' close(id) RDPSF = -1 RETURN END IF DO IPSTYP=1,MAXTYP I = NPARAM(IPSTYP, 1., CHECK, PAR, MAXPAR) IF ((LABEL .EQ. CHECK) .AND. (I .EQ. NPAR)) GO TO 1100 END DO write(*,*)'Inappropriate PSF: '//LABEL close(id) RDPSF = -1 RETURN C 1100 READ (id,301,IOSTAT=ISTAT) (PAR(I), I=1,NPAR) 301 FORMAT (1X, 6E13.6) IF (ISTAT .NE. 0) THEN write(*,*)'Error reading PSF.' RDPSF = -1 close(id) RETURN END IF NTERM = NEXP+NFRAC IF (NTERM .GE. 1) THEN DO K=1,NTERM READ(id,311,IOSTAT=ISTAT)((PSF(I,J,K),I=1,NPSF),J=1,NPSF) 311 FORMAT (1X, 6E13.6) IF (ISTAT .NE. 0) THEN close(id) RDPSF = -1 RETURN END IF END DO END IF close(id) RDPSF = 0 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 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 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 DAORAN (IDUM) C C RAN2 from Numerical Recipes C INTEGER IR(97) DATA IFF /0/ DATA M/714025/, IA/1366/, IC/150889/, RM/1.400511E-6/ 1000 CONTINUE IF ((IDUM .LT. 0) .OR. (IFF .EQ. 0)) THEN IFF = 1 IDUM = MOD(IABS(IC-IDUM), M) DO J=1,97 IDUM = MOD(IA*IDUM+IC,M) IR(J) = IDUM END DO IDUM = MOD(IA*IDUM+IC,M) IY = IDUM END IF J = 1+(97*IY)/M IF ((J .GT. 97) .OR. (J .LT. 1)) PAUSE IY = IR(J) DAORAN = IY*RM IDUM = MOD(IA*IDUM+IC,M) IR(J) = IDUM IF (DAORAN .LE. 0.) GO TO 1000 ! Stetson's modification C RETURN END! C C###################################################################### C REAL FUNCTION NRML (RANNUM) C C Convert a uniform probability distribution to a Gaussian distribution C with mean zero and standard deviation unity. C IMPLICIT NONE C REAL SQRT, ALOG C REAL P, RANNUM, SIGN, T 1000 P=RANNUM SIGN=-1. IF (P .GT. 0.5) THEN P=P-0.5 SIGN=1. ELSE IF (P .LE. 0.) THEN NRML = -1.E20 RETURN END IF T=SQRT(ALOG(1/P**2)) NRML=SIGN*(T- (2.30753+.27061*T) / (1.+T*(.99229+T*.04481)) ) RETURN 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####################################################################### 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.1, F8.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************************************* SUBROUTINE STRIP (ID, X, Y, MAG, SKY, SKIP, MAXSTR, . NSTAR, NDISAP, RADIUS, INDEX, HOLD) REAL X(MAXSTR), Y(MAXSTR), MAG(MAXSTR), SKY(MAXSTR) REAL HOLD(MAXSTR) INTEGER ID(MAXSTR), INDEX(MAXSTR) LOGICAL SKIP(MAXSTR) C NDISAP = 0 RADSQ = RADIUS**2 IF (NSTAR .LE. 1) RETURN DO I=1,NSTAR SKIP(I) = .FALSE. END DO C CALL QUICK (Y, NSTAR, INDEX) CALL RECTFYI (ID, NSTAR, INDEX, HOLD) CALL RECTFY (X, NSTAR, INDEX, HOLD) CALL RECTFY (MAG, NSTAR, INDEX, HOLD) CALL RECTFY (SKY, NSTAR, INDEX, HOLD) C DO 1200 I=1,NSTAR-1 IF (SKIP(I)) GO TO 1200 DO 1100 J=I+1,NSTAR IF (SKIP(J)) GO TO 1100 DY = Y(J)-Y(I) IF (DY .GT. RADIUS) GO TO 1200 DX = X(J)-X(I) IF (ABS(DX) .GT. RADIUS) GO TO 1100 IF (DX**2+DY**2 .GT. RADSQ) GO TO 1100 IF (MAG(J) .LE. MAG(I)) THEN SKIP(J) = .TRUE. GO TO 1100 ELSE SKIP(I) = .TRUE. GO TO 1200 END IF 1100 CONTINUE 1200 CONTINUE C ISTAR = 0 2000 CONTINUE IF (SKIP(NSTAR)) THEN NSTAR = NSTAR-1 NDISAP = NDISAP+1 GO TO 2000 END IF 2100 ISTAR = ISTAR+1 IF (ISTAR .GE. NSTAR) RETURN IF (SKIP(ISTAR)) THEN ID(ISTAR) = ID(NSTAR) X(ISTAR) = X(NSTAR) Y(ISTAR) = Y(NSTAR) MAG(ISTAR) = MAG(NSTAR) SKY(ISTAR) = SKY(NSTAR) SKIP(ISTAR) = .FALSE. NSTAR = NSTAR-1 NDISAP = NDISAP+1 GO TO 2000 ELSE GO TO 2100 END IF END C C####################################################################### C SUBROUTINE REGRP (ID, X, Y, MAG, SKY, CHI, DXOLD, DYOLD, . XCLAMP, YCLAMP, MAXSTR, NSTAR, FITRAD, LAST, INDEX, HOLD) C C======================================================================= C C This subroutine accepts a list of stellar coordinates, and C associates the stars into natural groups based on a critical C separation: stars within one critical separation of each other are C put into the same group; no star is within one critical separation of C any star outside its group. C C OFFICIAL DAO VERSION: 1985 August 15 C C====================================================================== C REAL X(MAXSTR), Y(MAXSTR), MAG(MAXSTR), SKY(MAXSTR) REAL CHI(MAXSTR), DXOLD(MAXSTR), DYOLD(MAXSTR) REAL XCLAMP(MAXSTR), YCLAMP(MAXSTR), HOLD(MAXSTR) INTEGER ID(MAXSTR), INDEX(MAXSTR) LOGICAL LAST(MAXSTR) C C----------------------------------------------------------------------- C C SECTION 1 C C Get set up. C IF (NSTAR .LE. 1) RETURN CRIT=2.*FITRAD CRITSQ=CRIT**2 C C Check that the stars are sorted by y-coordinate on input. C INDEX(1) = 1 DO I=2,NSTAR IF (Y(I) .GE. Y(I-1)) THEN INDEX(I) = I ELSE CALL QUICK (Y, NSTAR, INDEX) CALL RECTFY (X, NSTAR, INDEX, HOLD) GO TO 900 END IF END DO C 900 ITEST=0 ITOP=2 C C The stars are currently in a stack NSTAR stars long. The variable C ITEST will point to the position in the stack occupied by the star C which is currently the center of a circle of the critical radius, C within which we are looking for other stars; this also starts out C with a value of 1. ITOP points to the top position in the stack of C the stars which have not yet been assigned to groups; this starts C out with the value 2. Each time through, the C program goes down through the stack from ITOP and looks for stars C within the critical distance from the star at stack position C ITEST. When such a star is found, it changes places in the stack C with the star at ITOP and ITOP is incremented by one. When the C search reaches a star J such that Y(J)-Y(ITEST) > CRIT it is known C that no further stars will be found within the critical distance, the C pointer ITEST is incremented by one, and the search proceeds again C from the new value of ITOP. If the pointer ITEST catches up C with the pointer ITOP, that means that the group currently being C built up is complete. Then a new group is started beginning with C the star at the current position ITEST, ( = the instantaneous C value of ITOP), ITOP is incremented by 1, and the next group is built C up as before. C 2100 ITEST=ITEST+1 LAST(ITEST)=.FALSE. IF (ITEST .EQ. ITOP) THEN C C ITEST has reached ITOP; no other unassigned stars are within a C critical separation of any member of the current group. The group is C therefore complete. Signify this by setting LAST(ITEST-1)=.TRUE. C (ITEST = the current value of ITOP), and then increment the value of C ITOP by one. C J=ITEST-1 IF (J .GT. 0) LAST(J)=.TRUE. ITOP=ITOP+1 ! Increment ITOP C C If ITOP is greater than NSTAR at this point, then we are finished C (the last group contains one star). Otherwise, on with the search. C IF (ITOP .GT. NSTAR) THEN LAST(ITEST)=.TRUE. GO TO 3000 END IF END IF C C Now go through the list of unassigned stars, occupying positions ITOP C through NSTAR in the stack, to look for stars within the critical C distance of the star at position ITEST in the stack. If one is found, C move it up to stack position ITOP and increment ITOP by one. C XTEST=X(ITEST) YTEST=Y(ITEST) J=ITOP DO 2120 I=J,NSTAR DY=Y(I)-YTEST IF (DY .GT. CRIT) GO TO 2100 DX=X(I)-XTEST IF (ABS(DX) .GT. CRIT) GO TO 2120 IF (DX**2+DY**2 .GT. CRITSQ) GO TO 2120 C C This star is within the critical distance of the star at stack C position ITEST. Therefore it should be added to the current group by C moving it up to position ITOP in the stack, where the pointer ITEST C may eventually reach it. C CALL ASWAP (MAXSTR, ITOP, I, X, Y, INDEX) C C Now increment ITOP by 1 to point at the topmost unassigned star in the C stack. C ITOP=ITOP+1 C C If ITOP is greater than NSTAR, then all stars have been assigned to C groups, and we are finished. C IF (ITOP .GT. NSTAR) THEN DO K=ITEST,NSTAR-1 LAST(K)=.FALSE. END DO LAST(NSTAR)=.TRUE. GO TO 3000 END IF 2120 CONTINUE GO TO 2100 C 3000 CONTINUE C C Rectify the remaining quantities. C CALL RECTFYI (ID, NSTAR, INDEX, HOLD) CALL RECTFY (MAG, NSTAR, INDEX, HOLD) CALL RECTFY (SKY, NSTAR, INDEX, HOLD) CALL RECTFY (CHI, NSTAR, INDEX, HOLD) CALL RECTFY (DXOLD, NSTAR, INDEX, HOLD) CALL RECTFY (DYOLD, NSTAR, INDEX, HOLD) CALL RECTFY (XCLAMP, NSTAR, INDEX, HOLD) CALL RECTFY (YCLAMP, NSTAR, INDEX, HOLD) C RETURN END C C####################################################################### C SUBROUTINE ASWAP (MAXSTR, I, J, X, Y, INDEX) C C======================================================================= C C Make the I-th and J-th stars in the stack change places (J > I), C without otherwise altering the order of the stars. The other C arguments are self-evident. C C======================================================================= C IMPLICIT NONE INTEGER MAXSTR C REAL X(MAXSTR), Y(MAXSTR) INTEGER INDEX(MAXSTR) C REAL XHOLD, YHOLD INTEGER I, J, K, L, IHOLD C C----------------------------------------------------------------------- C c call ovrwrt ('ASWAP', 2) XHOLD=X(J) YHOLD=Y(J) IHOLD=INDEX(J) DO K=J,I+1,-1 L=K-1 X(K)=X(L) Y(K)=Y(L) INDEX(K)=INDEX(L) END DO X(I)=XHOLD Y(I)=YHOLD INDEX(I)=IHOLD RETURN END function indexpos(head,f1,n) character*80 head(*),f1*8 do 10 indexpos=1,n 10 if(head(indexpos)(1:8).eq.f1)return end