;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; This is the cfht_dph2psf source code in IDL. ; This program does the same thing as the C program cfht_dph2psf.c, that ; is it allows to compute PUEO's PSF associated with an object ; from the various dph (structure function) files saved during the ; acquisition and associated image of point sources. ; cfht_dph2psf take the name of a single parameter file as input, ; for example: cfht_dph2psf,'myparfile.par' ; The syntax of this parameter file must be the following: ; * PSFILE ; xxxxxxx <- Name of the output file for the calculated PSF ; ; * OBJECT ; xxxxxxx <- Name of the first dph file for the object ; ...... <- Name of other dph files for the object (optional) ; ; * REFERENCE 1 ; xxxxxxx <- Name of the image of the first reference ; xxxxxxx <- Name of the first dph file for the first reference ; ....... <- Name of other dph files for the first reference (optional) ; ; ....... ; ; * REFERENCE n ; xxxxxxx <- Name of the image of the nth reference ; xxxxxxx <- Name of the first dph file for the nth reference ; ....... <- Name of other dph files for the nth reference (optional) ; ; For more information, refer to the CFHT DPH2PSF manual which can be ; obtained from any of the following URL: ; http://www.cfht.hawaii.edu/~veran/cfht_dph2psf.html ; http://www.lyot.obspm.fr/~veran/cfht_dph2psf.html ; http://www.enst.fr/~veran/cfht_dph2psf.html ; ; -------------------------------------------------------------------------- ; ; Version 1.0 11 March 1997 Jean-Pierre Veran (veran@hplyot.obspm.fr) ; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; This function is the equivalent of the basename Unix facility. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; function basename,name return,strmid(name,rstrpos(name,'/')+1,1000) end ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; Return the content of the current line in parameter file ; If EOF is reached, return empty string ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; function get_line,line,fich line = '' if eof(fich) then return,line ; Return '' if eof readf,fich,line return,line end ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; Return next keyword value in parameter file. ; A keyword is a string starting with "* ". ; Return empty string if EOF is reached before finding any keyword. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; function find_next_keyword,line,fich while (not eof(fich)) do begin line = get_line(line,fich) line = strmid(line,0,8) ; limit keyword scope to 8 characters if strmid(line,0,2) eq '* ' then return,line endwhile return,line ; line has been set to empty string by get_line end ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; Read the FITS format file filename and put image data in the mono- ; dimensional floating point array a. ; Also check that the dimension are X_SIZE by Y_SIZE. ; Also if check is set and lambda is not zero, then check that the ; LAMBDA in the header is equal to lambda. Same for camsamp and CAMSAMP. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; pro read_fits_file,filename,a,lambda=lambda,camsamp=camsamp,check=check ; print,'Now reading ',filename a = readfits(filename,h,/silent) ;; read the NAXIS1 and NAXIS2 keyword to get image size naxes = fxpar(h,'NAXIS*',filename); if (naxes(0) ne 128) or (naxes(1) ne 128) then begin message,'ERROR: file '+filename+' is not a 128x128 image.' endif if keyword_set(check) then begin ;; read the LAMBDA and CAMSAMP keyword and check gotlambda = fxpar(h,'LAMBDA',filename) gotcamsamp = fxpar(h,'CAMSAMP',filename) if (lambda ne 0) then begin if gotlambda ne lambda then begin message,$ 'ERROR: The LAMBDA keyword must be the same for all the dph files' endif if gotcamsamp ne camsamp then begin message,$ 'ERROR: The CAMSAMP keyword must be the same for all the dph files' endif endif lambda = gotlambda; camsamp = gotcamsamp; endif end ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; This routine actually computes the PSF knowing the atmospheric structure ; function for the reference dphref, an image of the reference ref and ; the atmospheric structure function of the object. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; pro aob_dph2psfcomp,psf,dphobj,dphref,ref ; print,'Now computing PSF' ; Compute OTF of filter from difference of structure functions otffilt = exp(-0.5*(dphobj - dphref)) ; Mask OTF: All values at 0 in dph must be 0 in OTF except at origin w = where(dphobj eq 0) otffilt(w) = 0 otffilt(0) = 1 ; Compute OTF of reference image complexarr = fft(ref,-1) ; Multiply by OTF of filter complexarr = complexarr * otffilt ; Compute PSF psf = float(fft(complexarr,1)) ; Normalize psf = psf/total(psf) end ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; Main program. ; Returns 0 upon successful completion. ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; pro cfht_dph2psf,filename on_error,2 if n_params() ne 1 then begin printf,-1,'Usage: cfht_dph2psf parameter_file' return endif openr,fich,filename,/get_lun ; First keyword should be * PSFILE, giving the name of the output file if find_next_keyword(line,fich) ne '* PSFFIL' then begin message,'ERROR: Keyword * PSFFILE not found in '+filename endif if get_line(outfile,fich) eq '' then begin message,'ERROR: name of PSF file is not found' endif ; Create a basic header for result file psf = fltarr(128,128) mkhdr,hdr,float(psf) lambda = 0 & camsamp = 0 ; They will be read from FITS header ; Second keyword should be * OBJECT, giving the name of the dph for the object. if find_next_keyword(line,fich)ne '* OBJECT' then begin message,'ERROR: Keyword * OBJECT not found in '+filename endif dphobj = fltarr(128,128) n = 0 ; while (get_line(line,fich) ne '') do begin read_fits_file,line,tmp,lambda=lambda,camsamp=camsamp,/check dphobj=dphobj+tmp if n eq 0 then begin ;; Write the LAMBDA and CAMSAMP cards fxaddpar,hdr,'LAMBDA',lambda,'Central imaging wavelength' fxaddpar,hdr,'CAMSAMP',camsamp,'Size of one pixel (in arcsec)' endif ;; Write name of current dph into FITS header sn = strtrim(string(n+1),2) fxaddpar,hdr,'DPHOBJ'+sn,basename(line),$ 'Structure function '+sn+' for object' ;; Iterate n=n+1 endwhile if n eq 0 then begin message,'ERROR: at least one dph file must be given for the object' endif if (n gt 1) then dphobj=dphobj/n ;; Now we parse the description of the reference nref = 0 while(find_next_keyword(line,fich) eq "* REFERE") do begin sref = strtrim(string(nref+1),2) if get_line(line,fich) eq '' then begin message,'ERROR: File name for reference image '+sref+' is not found' endif read_fits_file,line,ref ;; Write name of current reference image in the header fxaddpar,hdr,'REF'+sref,basename(line),'Reference source number '+sref ;; Now look for the name of the corresponding dph files dphref=fltarr(128,128) n = 0 while (get_line(line,fich) ne '') do begin sn = strtrim(string(n+1),2) read_fits_file,line,tmp,lambda=lambda,camsamp=camsamp,/check dphref=dphref+tmp ;; Write name of current reference dph file in the header fxaddpar,hdr,'DREF'+sref+'_'+sn,basename(line),$ 'Structure function '+sn+' for reference '+sref ;; Iterate: next dph if any n=n+1 endwhile if n eq 0 then begin message,'ERROR: at least one dph file must be given for each dph' endif if n gt 1 then dphref=dphref/n aob_dph2psfcomp,tmp,dphobj,dphref,ref psf=psf+tmp ;; Iterate: next reference if any nref = nref + 1 endwhile if nref eq 0 then begin message,'ERROR: at least one reference must be used' endif if nref gt 1 then psf=psf/nref ;; write the PSF array to the FITS file writefits,outfile,psf,hdr ; print,psf(0,0),psf(1,0),psf(3,100),psf(64,64) end