;---------------------------------------------------------------------------------------------- ;---------------------------------------------------------------------------------------------- ; This code calculates the fractional transmission based on an ecam image of a standard star ; ;---------------------------------------------------------------------------------------------- ;---------------------------------------------------------------------------------------------- ;---------------------------------------------------------------------------------------------- ; Define Gaussian to be fit ;---------------------------------------------------------------------------------------------- function gfunct, x, y, P ;define Gaussian to be fit F=P[0]*exp(-((x-P[1])^2/(2*(P[3]^2))+(y-P[2])^2/(2*(P[3]^2)))) return, F end ;---------------------------------------------------------------------------------------------- ;---------------------------------------------------------------------------------------------- ; ----------------------------------------IMAGE SET 1------------------------------------------ ;---------------------------------------------------------------------------------------------- ;---------------------------------------------------------------------------------------------- ;---------------------------------------------------------------------------------------------- ; Define an image where the star was NOT in the slit, so the max flux ; can be determined and compared to when it is in the slit (the loop ; further down). ;---------------------------------------------------------------------------------------------- fraction=fltarr(42) ;number of images to be considered ****MAKE SURE TO UPDATE**** fracerr=fltarr(42) image=readfits('proc-e1267.fits') ;image of star not in slit ;image=readfits('proc-e0018.fits') image=image+32768.0 for i=0, size(image, /N_elements)-1 do begin if image[i] eq 0.0 then image[i]=1.0 endfor imageerr=sqrt(image) ;---------------------------------------------------------------------------------------------- ; Find max and min pixel values of image ; Fit guassian function ; b = the x position of the max value pixel. The data is read as a ; list of values, so to find the row you need to divide by 511 (the ; number of pixels in the row) ; c = same as b except for the y position ; aper, image, c, b, flux, errflux, sky, skyerr, 1, [30], [50, 100], [0.0,30000.0], /flux ; aper = aperture function defined in idl. It will draw an aperature ; of a defined size. Here it's defined as 30 pixels in radius (it ; assumes a circular aperture by default) ; ; flux = an idl routine that calculate the integrated flux of an ; aperture that was define ; ; sky = idl routine that find the integrated flux of a sky area same ; aperatue size as the previous aperature ; ; 0.0, 30000 sets the min and max value it will consider. This is a ; problem for stars that are faint because the cosmic rays will be ; brighter and not necessarily 30,000. You can set this down based on ; the higherst point of your star. ; ; [50, 100] is the annulus of the sky. The ring will go around the ; star from 50 though 100 pixels out from the star's central location. ; ; In order to try and sort out the cosmic rays from the star, try ; making a list of the top 15 max values or so and their location. ; Then check where the pixel values are in relation to each other. If ; it's a cosmic ray, it will only be 1-4, but a star will have them ; all together. Or you can do an aperature and integrated flux for ; each of the top 15 max positions and choose the one with the highest ; integrated flux to be your star. Just make sure to print out the ; max position so you can double check in the image that the spot is accurate. ;---------------------------------------------------------------------------------------------- xr=findgen(511) ;# define pixels in x direction for image yc=findgen(511) ;# define pixels in y direction for image ;------------------------------------------------------------------------------------------------ ; Defining position of star manually if it can't be found with a search. Find peak pixel in DS9 ;------------------------------------------------------------------------------------------------ x = 35 y = 297 aper, image, x, y, flux, errflux, sky, skyerr, 1, [8], [50, 100], [0.0,3000.0], /flux ;b=floor(where(image eq max(image))/511.0) ;c=(where(image eq max(image))/511.0-floor(where(image eq max(image))/511.0))*511.0 parinfo = replicate({value:0.D, fixed:0, limited:[0,0], limits:[0.0,0.0]},4) parinfo[1].fixed=1 parinfo[2].fixed=1 parinfo[*].value=[8000.0,x,y,15.0] P=mpfit2dfun('gfunct', xr, yc, image, imageerr, PARINFO=parinfo, NITER=niter, PERROR=perror, STATUS=status, BESTNORM=bestnorm, QUIET=0) print, 'P=',P fit=fltarr(size(yc, /N_elements)) fit=8000.0*exp(-((yc-224.0)^2/(2*(6.8^2)))) ;---------------------------------------------------------------------------------------------- ; Apply to other ecam images in the same set. ; Define the first ecam number for the set of ecam images. I.e., if ; the first figure is called proc-e0042.fits and the last is ; proc-e0097.fits, then start_img = 42 and end_img = 97 ; The total number of images being considered is than 97 - 42 + 1 ;---------------------------------------------------------------------------------------------- start_img = 1281 end_img = 1322 tot_img = end_img - start_img + 1 for i=start_img, end_img do begin print, i image2=readfits('proc-e'+strn(i)+'.fits') image2=image2+32768.0 x = 201 ;manually defining center of star in x and y y = 218 ; b=floor(where(image2 eq max(image2))/511.0) ; c=(where(image2 eq max(image2))/511.0-floor(where(image2 eq max(image2))/511.0))*511.0 aper, image2, x, y, slitflux, slitfluxerr, skyslit, skysliterr, 1, [16], [50, 100], [0.0,3000.0], /flux fraction[i-start_img]=slitflux/flux fracerr[i-start_img]=fraction*sqrt((slitfluxerr/slitflux)^2.0+(errflux/flux)^2.0) print, fraction[i-start_img] endfor fraction=fraction(where(finite(fraction) eq 1)) print, fraction med=median(fraction) med=1-med ;actually want 1-transmission to get light in slit average=total(fraction)/size(fraction, /N_elements) average=1-average averagevector=replicate(average, tot_img) aveerrsquared=0 for i=0, size(fraction, /N_elements)-1 do begin aveerrsquared=aveerrsquared+fracerr[i]^2.0 endfor aveerr=sqrt(aveerrsquared)/size(fraction, /N_elements) sigma=stddev(fraction) sigmavector=replicate(sigma, tot_img) plot, 1-fraction, psym=1 oplot, averagevector-sigmavector, linestyle=5 oplot, averagevector+sigmavector, linestyle=5 print, 'median=', med ;Given the movement of the star, median will be more accurate than average print, 'average=', average, '+/-', aveerr print, 'stddev=', sigma print, 'percent error=', sigma/average close,1 openw, 1, 'fluxcal1.txt' for i=0, size(fraction, /N_elements)-1 do begin printf, 1, i, string(9b), 1-fraction[i], string(9b), averagevector[i], string(9b), sigmavector[i] endfor close,1 ;---------------------------------------------------------------------------------------------- ; Values are calculated based on a set of ecam images that were taken ; for the exposure of the star. If you want another set of ecam ; images to be analyzed for another stellar observations, type 1 to ; continue, else the program will end. ;---------------------------------------------------------------------------------------------- print, 'enter 1 to continue' read, continue if continue eq 1 then begin endif else begin stop endelse ;---------------------------------------------------------------------------------------------- ;---------------------------------------------------------------------------------------------- ; ----------------------------------------IMAGE SET 2------------------------------------------ ;---------------------------------------------------------------------------------------------- ;---------------------------------------------------------------------------------------------- ;---------------------------------------------------------------------------------------------- ;---------------------------------------------------------------------------------------------- ; Start again with a new set of ecam images and go through all steps ;---------------------------------------------------------------------------------------------- ;---------------------------------------------------------------------------------------------- fraction=fltarr(23) image=readfits('proc-e0149.fits') image=image+32768.0 xr=findgen(511) yc=findgen(511) ;b=floor(where(image eq max(image))/511.0) ;c=(where(image eq max(image))/511.0-floor(where(image eq max(image))/511.0))*511.0 x = 33 y = 300 ;x = 41 ;y = 291 aper, image, x, y, flux, errflux, sky, skyerr, 1, [20], [50, 100], [0.0,3000.0], /flux plot, yc, image[*,y], xrange=[0,100] ;--------------------------------------------------------------------------------------- ; List of ecam images with star in slit ;-------------------------------------------------------------------------------------- start_img =88 end_img = 110 tot_img = end_img - start_img + 1 for i=start_img, end_img do begin print, i stringnumber = strn(i,FORMAT='(I04)') ; makes the image number an integer (I) of 4 characters, ; with all leading characters set to 0. If i = 88, then ; stringnumber is 0088, and if i = 100, then stringnumber = 0100. image2 = readfits('proc-e'+stringnumber+'.fits') image2=image2+32768.0 ; if i lt 100 then image2=readfits('proc-e00'+strn(i)+'.fits') ; if i gt 99 then image2=readfits('proc-e0'+strn(i)+'.fits') ; This way also works to deal with the 0 issue ;b=floor(where(image2 eq max(image2))/511.0) ;c=(where(image2 eq max(image2))/511.0-floor(where(image2 eq max(image2))/511.0))*511.0 x = 199 y = 216 aper, image2, x, y, slitflux, slitfluxerr, skyslit, skysliterr, 1, [20], [50, 100], [0.0,3000.0], /flux fraction[i-start_img]=slitflux/flux print, fraction[i-start_img] endfor fraction=fraction(where(finite(fraction) eq 1)) print, fraction med=median(fraction) med=1-med ;actually want 1-transmission to get light in slit average=total(fraction)/size(fraction, /N_elements) average=1-average averagevector=replicate(average, tot_img) sigma=stddev(fraction) sigmavector=replicate(sigma, tot_img) plot, 1-fraction, psym=1 oplot, averagevector-sigmavector, linestyle=5 oplot, averagevector+sigmavector, linestyle=5 print, 'median=', med print, 'average=', average print, 'stddev=', sigma print, 'percent error=', sigma/average openw, 1, 'fluxcal2.txt' for i=0, size(fraction, /N_elements)-1 do begin printf, 1, i, string(9b), 1-fraction[i], string(9b), averagevector[i], string(9b), sigmavector[i] endfor close,1 print, 'enter 1 to continue' read, continue if continue eq 1 then begin endif else begin stop endelse ;---------------------------------------------------------------------------------------------- ;---------------------------------------------------------------------------------------------- ; ----------------------------------------IMAGE SET 3------------------------------------------ ;---------------------------------------------------------------------------------------------- ;---------------------------------------------------------------------------------------------- ;---------------------------------------------------------------------------------------------- ; Start again with a new set of ecam images and go through all steps ;---------------------------------------------------------------------------------------------- Print, "Image set 3, image proc-e0149.fits \n" image=readfits('proc-e0149.fits') image=image+32768.0 xr=findgen(511) yc=findgen(511) x = 33 y = 300 ;b=floor(where(image eq max(image))/511.0) ;c=(where(image eq max(image))/511.0-floor(where(image eq max(image))/511.0))*511.0 aper, image, x, y, flux, errflux, sky, skyerr, 1, [15], [50, 100], [0.0,3000.0], /flux plot, yc, image[*,y], xrange=[190,240] ;---------------------------------------------------------------------------------------------- ; Loop through images where star is on slit ;---------------------------------------------------------------------------------------------- READCOL,'position_A0V0003.dat',image_numb,x,y,FORMAT='A,I,I' ;Read 3 column data file with integer ;columns, assigned to the image number ;and the x and y position fraction=fltarr(N_ELEMENTS(image_numb)) start_img = 160 end_img = 203 tot_img = end_img - start_img + 1 for i=0, 43 do begin print, i print, image_numb[i] ; image_numb[i] = strn(image_numb[i],FORMAT='(I04)') ; makes the image number an integer (I) of 4 characters, ; with all leading characters set to 0. If i = 88, then ; stringnumber is 0088, and if i = 100, then stringnumber = 0100. image2 = readfits('proc-e'+image_numb[i]+'.fits') image2=image2+32768.0 aper, image2, x[i], y[i], slitflux, slitfluxerr, skyslit, skysliterr, 1, [15], [50, 100], [0.0,300.0], /flux fraction[i]=slitflux/flux print, fraction[i] endfor fraction=fraction(where(finite(fraction) eq 1)) print, fraction med=median(fraction) med=1-med ;actually want 1-transmission to get light in slit average=total(fraction)/size(fraction, /N_elements) average=1-average averagevector=replicate(average, tot_img) sigma=stddev(fraction) sigmavector=replicate(sigma, tot_img) plot, 1-fraction, psym=1 oplot, averagevector-sigmavector, linestyle=5 oplot, averagevector+sigmavector, linestyle=5 print, 'median=', med print, 'average=', average print, 'stddev=', sigma print, 'percent error=', sigma/average openw, 1, 'fluxcal3.txt' for i=0, size(fraction, /N_elements)-1 do begin printf, 1, i, string(9b), 1-fraction[i], string(9b), averagevector[i], string(9b), sigmavector[i] endfor close,1 print, 'enter 1 to continue' read, continue if continue eq 1 then begin endif else begin stop endelse ;---------------------------------------------------------------------------------------------- ;---------------------------------------------------------------------------------------------- ; ----------------------------------------IMAGE SET 4------------------------------------------ ;---------------------------------------------------------------------------------------------- ;---------------------------------------------------------------------------------------------- ;---------------------------------------------------------------------------------------------- ; Start again with a new set of ecam images and go through all steps ;---------------------------------------------------------------------------------------------- Print, "Image set 4, image proc-e0266.fits \n" image=readfits('proc-e0266.fits') image=image+32768.0 xr=findgen(511) yc=findgen(511) x = 90 y = 261 ;b=floor(where(image eq max(image))/511.0) ;c=(where(image eq max(image))/511.0-floor(where(image eq max(image))/511.0))*511.0 aper, image, x, y, flux, errflux, sky, skyerr, 1, [18], [50, 100], [0.0,30000.0], /flux plot, yc, image[*,y], xrange=[190,240] ;---------------------------------------------------------------------------------------------- ; Loop through images where star is on slit ;---------------------------------------------------------------------------------------------- READCOL,'position_A0V0004.dat',image_numb,x,y,FORMAT='A,I,I' ;Read 3 column data file with integer ;columns, assigned to the image number ;and the x and y position fraction=fltarr(N_ELEMENTS(image_numb)) start_img = 279 end_img = 300 tot_img = end_img - start_img + 1 for i=0, N_ELEMENTS(image_numb)-1 do begin print, i print, image_numb[i] ; image_numb[i] = strn(image_numb[i],FORMAT='(I04)') ; makes the image number an integer (I) of 4 characters, ; with all leading characters set to 0. If i = 88, then ; stringnumber is 0088, and if i = 100, then stringnumber = 0100. image2 = readfits('proc-e'+image_numb[i]+'.fits') image2=image2+32768.0 aper, image2, x[i], y[i], slitflux, slitfluxerr, skyslit, skysliterr, 1, [19], [50, 100], [0.0,30000.0], /flux fraction[i]=slitflux/flux print, fraction[i] endfor fraction=fraction(where(finite(fraction) eq 1)) print, fraction med=median(fraction) med=1-med ;actually want 1-transmission to get light in slit average=total(fraction)/size(fraction, /N_elements) average=1-average averagevector=replicate(average, tot_img) sigma=stddev(fraction) sigmavector=replicate(sigma, tot_img) plot, 1-fraction, psym=1 oplot, averagevector-sigmavector, linestyle=5 oplot, averagevector+sigmavector, linestyle=5 print, 'median=', med print, 'average=', average print, 'stddev=', sigma print, 'percent error=', sigma/average openw, 1, 'fluxcal4.txt' for i=0, size(fraction, /N_elements)-1 do begin printf, 1, i, string(9b), 1-fraction[i], string(9b), averagevector[i], string(9b), sigmavector[i] endfor close,1 print, 'enter 1 to continue' read, continue if continue eq 1 then begin endif else begin stop endelse stop ; read through images one at a time ; read in a data file with the x and y value rather then looking for ; them or specifying them. start_img = 160 end_img = 201 tot_img = end_img - start_img + 1 fraction=fltarr(tot_img) for i=start_img, end_img do begin print, i stringnumber = strn(i,FORMAT='(I04)') ; makes the image number an integer (I) of 4 characters, ; with all leading characters set to 0. If i = 88, then ; stringnumber is 0088, and if i = 100, then stringnumber = 0100. image2 = readfits('proc-e'+stringnumber+'.fits') image2=image2+32768.0 ; if i lt 100 then image2=readfits('proc-e00'+strn(i)+'.fits') ; if i gt 99 then image2=readfits('proc-e0'+strn(i)+'.fits') ; This way also works to deal with the 0 issue ;b=floor(where(image2 eq max(image2))/511.0) ;c=(where(image2 eq max(image2))/511.0-floor(where(image2 eq max(image2))/511.0))*511.0 x = 201 y = 217 aper, image2, x, y, slitflux, slitfluxerr, skyslit, skysliterr, 1, [19], [50, 100], [0.0,3000.0], /flux fraction[i-start_img]=slitflux/flux print, fraction[i-start_img] endfor fraction=fraction(where(finite(fraction) eq 1)) print, fraction med=median(fraction) med=1-med ;actually want 1-transmission to get light in slit average=total(fraction)/size(fraction, /N_elements) average=1-average averagevector=replicate(average, tot_img) sigma=stddev(fraction) sigmavector=replicate(sigma, tot_img) plot, 1-fraction, psym=1 oplot, averagevector-sigmavector, linestyle=5 oplot, averagevector+sigmavector, linestyle=5 print, 'median=', med print, 'average=', average print, 'stddev=', sigma print, 'percent error=', sigma/average openw, 1, 'fluxcal3.txt' for i=0, size(fraction, /N_elements)-1 do begin printf, 1, i, string(9b), 1-fraction[i], string(9b), averagevector[i], string(9b), sigmavector[i] endfor close,1 print, 'enter 1 to continue' read, continue if continue eq 1 then begin endif else begin stop endelse ;HILT600_4 ;---------------------------------------------------------------------------------------------- ; Start again with a new set of ecam images and go through all steps ;---------------------------------------------------------------------------------------------- fraction=fltarr(56) image=readfits('proc-e1008.fits') image=image+32768.0 xr=findgen(511) yc=findgen(511) b=floor(where(image eq max(image))/511.0) c=(where(image eq max(image))/511.0-floor(where(image eq max(image))/511.0))*511.0 aper, image, c, b, flux, errflux, sky, skyerr, 1, [30], [50, 100], [0.0, 300.0], /flux plot, yc, image[*,b], xrange=[230,330] for i=1009, 1064 do begin print, i image2=readfits('proc-e'+strn(i)+'.fits') image2=image2+32768.0 b=floor(where(image2 eq max(image2))/511.0) c=(where(image2 eq max(image2))/511.0-floor(where(image2 eq max(image2))/511.0))*511.0 aper, image2, c, b, slitflux, slitfluxerr, skyslit, skysliterr, 1, [30], [50, 100], [0.0, 300.0], /flux fraction[i-1009]=slitflux/flux print, fraction[i-1009] endfor fraction=fraction(where(finite(fraction) eq 1)) print, fraction med=median(fraction) med=1-med ;actually want 1-transmission to get light in slit average=total(fraction)/size(fraction, /N_elements) average=1-average averagevector=replicate(average, 56) sigma=stddev(fraction) sigmavector=replicate(sigma, 56) plot, 1-fraction, psym=1 oplot, averagevector-sigmavector, linestyle=5 oplot, averagevector+sigmavector, linestyle=5 print, 'median=', med print, 'average=', average print, 'stddev=', sigma print, 'percent error=', sigma/average openw, 1, 'fluxcal4.txt' for i=0, size(fraction, /N_elements)-1 do begin printf, 1, i, string(9b), 1-fraction[i], string(9b), averagevector[i], string(9b), sigmavector[i] endfor close,1 print, 'enter 1 to continue' read, continue if continue eq 1 then begin endif else begin stop endelse ;HILT600_5 ;---------------------------------------------------------------------------------------------- ; Start again with a new set of ecam images and go through all steps ;---------------------------------------------------------------------------------------------- fraction=fltarr(56) image=readfits('proc-e1268.fits') image=image+32768.0 xr=findgen(511) yc=findgen(511) b=floor(where(image eq max(image))/511.0) c=(where(image eq max(image))/511.0-floor(where(image eq max(image))/511.0))*511.0 aper, image, c, b, flux, errflux, sky, skyerr, 1, [30], [50, 100], [0.0,25000.0], /flux plot, yc, image[*,b], xrange=[230,330] for i=1269, 1324 do begin print, i image2=readfits('proc-e'+strn(i)+'.fits') image2=image2+32768.0 b=floor(where(image2 eq max(image2))/511.0) c=(where(image2 eq max(image2))/511.0-floor(where(image2 eq max(image2))/511.0))*511.0 aper, image2, c, b, slitflux, slitfluxerr, skyslit, skysliterr, 1, [30], [50, 100], [0.0,25000.0], /flux fraction[i-1269]=slitflux/flux print, fraction[i-1269] endfor fraction=fraction(where(finite(fraction) eq 1)) print, fraction med=median(fraction) med=1-med ;actually want 1-transmission to get light in slit average=total(fraction)/size(fraction, /N_elements) average=1-average averagevector=replicate(average, 57) sigma=stddev(fraction) sigmavector=replicate(sigma, 57) plot, 1-fraction, psym=1 oplot, averagevector-sigmavector, linestyle=5 oplot, averagevector+sigmavector, linestyle=5 print, 'median=', med print, 'average=', average print, 'stddev=', sigma print, 'percent error=', sigma/average openw, 1, 'fluxcal5.txt' for i=0, size(fraction, /N_elements)-1 do begin printf, 1, i, string(9b), 1-fraction[i], string(9b), averagevector[i], string(9b), sigmavector[i] endfor close,1 end ;------------------------------------------------------------------------------------------------ ; CUT STUFF ;xr=findgen(511) ;# pixels in x direction ;yc=findgen(511) ;# pixels in y direction ;list = SORT(image) ; sorts the values from max to min, but only prints their xy positions ;maximums = image[REVERSE(SORT(image))] ; sorts only the pixels value (don't think it retains position) ;list_short_x = fltarr(550) ;list_short_y = fltarr(550) ;max_short = fltarr(550) ;brightness = fltarr(550) ;for i=0, 549 do begin ; aper, image, x, y, flux, errflux, sky, skyerr, 1, [30], [50, 100], [0.0,30000.0], /flux ; list_short_x[i] = list[i] ; prints the pixel number (starting from 0 - 511*511) ; list_short_y[i] = list[i] - floor(list[i]) ; list_short_x[i] = floor(list[i]/511.0)+1 ; list_short_y[i] = (list[i]/511.0 - floor(list[i]/511.0))*511.0 ; max_short[i] = maximums[i] ; print, list_short_x[i], list_short_y[i], max_short[i] ; aper, image, list_short_x[i], list_short_y[i], flux, errflux, sky, skyerr, 1, [30], [50, 100], [0.0,30000.0], /flux ; brightness[i] = flux ;endfor ;-------------------------------------------------------------------------------------------- ;;plot, yc, image[*,b], xrange=[210,240] ;oplot, yc, P[0]*exp(-((yc-P[2])^2/(2*(P[3]^2)))) ;oplot, yc, fit, linestyle=5 ;---------------------------------------------------------------------------------------------- ; Apply to other ecam images in the same set. ; Define the first ecam number for the set of ecam images. I.e., if ; the first figure is called proc-e0042.fits and the last is ; proc-e0097.fits, then start_img = 42 and end_img = 97 ; The total number of images being considered is than 97 - 42 + 1 ;for i=42, 97 do begin ; print, i ; image2=readfits('proc-e00'+strn(i)+'.fits') ; image2=image2+32768.0 ; b=floor(where(image2 eq max(image2))/511.0) ; c=(where(image2 eq max(image2))/511.0-floor(where(image2 eq max(image2))/511.0))*511.0 ; aper, image2, c, b, slitflux, slitfluxerr, skyslit, skysliterr, 1, [30], [50, 100], [0.0,30000.0], /flux ; fraction[i-42]=slitflux/flux ; fracerr[i-42]=fraction*sqrt((slitfluxerr/slitflux)^2.0+(errflux/flux)^2.0) ; print, fraction[i-42] ;endfor ;fraction=fraction(where(finite(fraction) eq 1)) ;print, fraction ;med=median(fraction) ;med=1-med ;actually want 1-transmission to get light in slit ;average=total(fraction)/size(fraction, /N_elements) ;average=1-average ;averagevector=replicate(average, 56) ;aveerrsquared=0 ;for i=0, size(fraction, /N_elements)-1 do begin ;aveerrsquared=aveerrsquared+fracerr[i]^2.0 ;endfor ;aveerr=sqrt(aveerrsquared)/size(fraction, /N_elements) ;sigma=stddev(fraction) ;sigmavector=replicate(sigma, 56) ;----------------------------------------------------------------------------------------------