;+ ; fm.pro IDL procedure ; ; PURPOSE ; Compute absolute/excess coronal mass. ; Read FITS images and define a region of interest. ; ; SYNTAX ; fm, 'cme_fits_file' ; fm, 'cme_fits_file', 'ref_fits_file' ; ; PROCEDURES USED ; rfits ; fxpar ; defroi ; xmass ; ; HISTORY ; Author: Andrew L. Stanger 2 June 1998 ;- PRO fm, cme_fits_file, ref_fits_file ;--- Initialize variables. mx = 0 ; =0 cme image only, =1 cme & ref images tscope = 0 ; Telescope: =1, SMM C/P, =2, MK3. xdim_cme = 0 ; X-axis dimension ydim_cme = 0 ; Y-axis dimension nbypix_cme = 0 ; # bytes/pixel xcen_cme = 0 ; X-axis center ycen_cme = 0 ; Y-axis center crradius_cme = 0.0 ; # pixels/Rsun telescop_cme = '' ; Telescope xdim_ref = 0 ; X-axis dimension ydim_ref = 0 ; Y-axis dimension nbypix_ref = 0 ; # bytes/pixel xcen_ref = 0 ; X-axis center ycen_ref = 0 ; Y-axis center crradius_ref = 0.0 ; # pixels/Rsun telescop_ref = '' ; Telescope ;--- Read CME FITS image. imgcme = rfits (cme_fits_file, header=hdu_cme, /scale) dspcme = imgcme xdim_cme = fxpar (hdu_cme, "NAXIS1") ydim_cme = fxpar (hdu_cme, "NAXIS2") nbypix_cme = fxpar (hdu_cme, "BITPIX") * 8 xcen_cme = fxpar (hdu_cme, "CRPIX1") ycen_cme = fxpar (hdu_cme, "CRPIX2") crradius_cme = fxpar (hdu_cme, "CRRADIUS") pixrs = crradius_cme ; print, 'fr> crradius: ', crradius, ' pixrs: ', pixrs dispmin_cme = fxpar (hdu_cme, "DISPMIN") dispmax_cme = fxpar (hdu_cme, "DISPMAX") print, 'CME dispmin/max: ', dispmin_cme, dispmax_cme telescop_cme = fxpar (hdu_cme, "TELESCOP") IF (strpos (telescop_cme, "SMM C/P") NE -1) THEN $ ; SMM C/P image. BEGIN tscope = 1 dspcme = imgcme ; END IF (strpos (telescop_cme, "MK3") NE -1) THEN $ ; MK3 image. BEGIN tscope = 2 imgvig = fltarr (512, 512) openr, 1, '/home/stanger/idl/lib/icat/mk3vig.template' readu, 1, imgvig close, 1 dspcme = imgcme * imgvig ; Apply "vignetting" print, 'Apply MK3 vignetting to dspcme.' END ; date = fxpar (hdu_cme, "DATE-OBS") ; time = fxpar (hdu_cme, "TIME-OBS") ; date_time = date + ' ' + time ; solar_info = pb0r (date_time,/arcsec) ; Get solar P, B0, Rsun. ; Pangle = solar_info [0] ; P-angle ; B0 = solar_info [1] ; B0 angle ; Rsunarc = solar_info [2] ; Solar radius (arcsec). window, xsize=xdim_cme, ysize=ydim_cme ; Create image window. ;--- Reference Image ? IF (N_PARAMS() EQ 2) THEN $ ; Reference (pre-event) image supplied. BEGIN ;{ imgref = rfits (ref_fits_file, header=hdu_ref, /scale) dspref = imgref xdim_ref = fxpar (hdu_ref, "NAXIS1") ydim_ref = fxpar (hdu_ref, "NAXIS2") nbypix_ref = fxpar (hdu_ref, "BITPIX") * 8 xcen_ref = fxpar (hdu_ref, "CRPIX1") ycen_ref = fxpar (hdu_ref, "CRPIX2") crradius_ref = fxpar (hdu_ref, "CRRADIUS") telescop_ref = fxpar (hdu_ref, "TELESCOP") err = 0 IF (xdim_cme NE xdim_ref) THEN err = -1 IF (ydim_cme NE ydim_ref) THEN err = -2 IF (nbypix_cme NE nbypix_ref) THEN err = -3 IF (xcen_cme NE xcen_ref) THEN err = -4 IF (ycen_cme NE ycen_ref) THEN err = -5 IF (err NE 0) THEN RETURN if (strpos (telescop_ref, "SMM C/P") NE -1) THEN $ ; SMM C/P image. BEGIN tscope = 1 dspref = imgref END if (strpos (telescop_ref, "MK3") NE -1) THEN $ ; MK3 image. BEGIN tscope = 2 imgvig = fltarr (512, 512) openr, 1, '/home/stanger/idl/lib/icat/mk3vig.template' readu, 1, imgvig close, 1 dspref = imgref * imgvig ; Apply "vignetting" print, 'Apply MK3 vignetting to dspref.' END mx = 1 lct, '/home/stanger/color/dif.lut' ; Load difference color map. dspdif = dspcme - dspref imin = MIN (dspdif, max=imax) print, 'dif min/max: ', imin, imax ;imin = imin * 1.5 imin = imin * 1.5 ;imax = imax * 0.5 imax = imax * 0.005 ;adjust this number to brighten difference image if (abs (imin) LT abs (imax)) THEN imax = - imin if (abs (imin) GT abs (imax)) THEN imin = - imax ; print, 'dif min/max: ', imin, imax dspdif = bytscl (dspdif, min=imin, max=imax, top=250) tv, dspdif ; Display difference image. END $ ;} ;--- No Reference image. ELSE $ ; No reference image supplied. BEGIN ;{ mx = 0 imgref = 0 lct, '/home/stanger/color/bwy.lut' ; Load blue-white-yellow color map. imin = MIN (dspcme, max=imax) ; print, 'CME min/max range: ', imin, imax IF (imin LT dispmin_cme) THEN imin = dispmin_cme IF (imax GT dispmax_cme) THEN imax = dispmax_cme ; print, 'CME min/max range: ', imin, imax ; print, 'CME disp min/max : ', dispmin_cme, dispmax_cme IF (dispmin_cme NE dispmax_cme) THEN $ BEGIN dspcmeb = bytscl (dspcme, min=imin, max=imax, top=249) END $ ELSE $ dspcmeb = bytscl (dspcme) tv, dspcmeb ; Display image. ; tvscl, dspcme END ;} print, "Define a polygon area for mass integration using the mouse." ; print, "Click with left mouse button to define polygon vertices." ; print, "Click with right mouse button to terminate the polygon." ;--- Define Region-of-interest (ROI). ;RETURN polypix = DEFROI (xdim_cme, ydim_cme, xv, yv,/RESTORE) ; array of pixel coordinates. ;PLOTS,xv,yv,COLOR=255,/DEVICE PLOTS,xv,yv,COLOR=140,/DEVICE size_roi = SIZE (polypix) nloc = size_roi [1] xloc = intarr (nloc) yloc = intarr (nloc) ;--- Convert 1-D coordinates to 2-D coordinates. FOR i = 0, nloc-1 DO $ BEGIN yloc [i] = polypix (i) / xdim_cme xloc [i] = polypix (i) - yloc (i) * xdim_cme END size_xv = SIZE (xv) nv = size_xv [1] ; print, 'size [xv]: ', size (xv) ; print, 'size [yv]: ', size (yv) ; FOR i = 0, nv - 1 DO $ ; print, 'i, xv, yv: ', i, xv [i], yv [i] ; print, 'size (polypix): ', size (polypix) ; FOR i = 0, nloc - 1 DO $ ; print, 'polypix, xloc, yloc: ', polypix (i), xloc (i), yloc (i) ;--- Compute integrated absolute/excess mass inside polygon. xm = xmass (imgref, imgcme, xdim_cme, ydim_cme, nbypix, $ xcen_cme, ycen_cme, pixrs, xloc, yloc, nloc, tscope) ; ********************* LOOK HERE ANDY ********************************** ;--- Compute integrated mass inside reference box in image ref_nloc=10000 ;Number of pixels in a 100x100 region ref_xloc=intarr(ref_nloc) ;x-axis array ref_yloc=intarr(ref_nloc) ;y-axis array ref_xvalues=indgen(100)+200 ;look at pixels 200-299 (x-axis) ref_yvalues=indgen(100)+380 ;look at pixels 380-480 (y-axis) ;---Fill arrays for xmass.pro FOR i=0, (ref_nloc)-1, 100 DO BEGIN ref_xloc(i:i+99)=ref_xvalues ref_yloc(i:i+99)=ref_yvalues ENDFOR ref_xm = xmass (imgref, imgcme, xdim_cme, ydim_cme, nbypix, $ xcen_cme, ycen_cme, pixrs, ref_xloc, ref_yloc, ref_nloc, tscope) ;--- End of reference box computation ; ************************** END **************************************** ;xm=xm * (10/3) ; Correction factor brought up by Joan and Art print, 'xm: ', xm print, 'ref_xm: ',ref_xm ;new_img=it_annotate(tvrd(),hdu_cme,hdu_ref) RETURN END