;+ ; fxm.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 fxm, cme_fits_file, ref_fits_file ;--- Initialize variables. i = 0L nloc = 0L red = bytarr (256) green = bytarr (256) blue = bytarr (256) 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") cmedate = fxpar (hdu_cme, "DATE-OBS") cmetime = fxpar (hdu_cme, "TIME-OBS") print, 'cmedate: ', cmedate print, 'cmetime: ', cmetime cmeday = strmid (cmedate, 0, 2) cmemonth = strmid (cmedate, 3, 2) cmeyear = strmid (cmedate, 6, 2) cmehour = strmid (cmetime, 0, 2) cmeminute = strmid (cmetime, 3, 2) cmesecond = strmid (cmetime, 6, 2) print, 'cme year, month, day: ', cmeyear, cmemonth, cmeday print, 'cme hour, minute, second: ', cmehour, cmeminute, cmesecond 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/vignet/mk3vig.template' readu, 1, imgvig close, 1 dspcme = imgcme * imgvig ; Apply "vignetting" print, 'Apply MK3 vignetting to dspcme.' END ; 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") refdate = fxpar (hdu_ref, "DATE-OBS") reftime = fxpar (hdu_ref, "TIME-OBS") refday = strmid (refdate, 0, 2) refmonth = strmid (refdate, 3, 2) refyear = strmid (refdate, 6, 2) refhour = strmid (reftime, 0, 2) refminute = strmid (reftime, 3, 2) refsecond = strmid (reftime, 6, 2) aoi_file = cmeyear + cmemonth + cmeday + '.' + $ cmehour + cmeminute + '-' + $ refhour + refminute + '.aoi.gif' log_file = cmeyear + cmemonth + cmeday + '.' + $ cmehour + cmeminute + '-' + $ refhour + refminute + '.fm' 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/vignet/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 ; IF (tscope EQ 1) THEN imax = imax * 0.5 IF (tscope EQ 2) THEN imax = imax * 0.02 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 aoi_file = cmeyear + cmemonth + cmeday + '.' + $ cmehour + cmeminute + '.aoi.gif' log_file = cmeyear + cmemonth + cmeday + '.' + $ cmehour + cmeminute + '.fm' END ;} print, "Create a polygon for mass integration by clicking with 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 an array of pixel coordinates. ;RETURN polypix = DEFROI (xdim_cme, ydim_cme, xv, yv, /NOFILL) size_roi = SIZE (polypix) nloc = size_roi [1] xloc = intarr (nloc) yloc = intarr (nloc) ; print, 'size (polypix): ', size_roi ; print, 'nloc: ', nloc ;--- Convert 1-D coordinates to 2-D coordinates. FOR i = 0L, 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 = 0L, 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 (log_file, ref_fits_file, cme_fits_file, imgref, imgcme, $ xdim_cme, ydim_cme, nbypix, xcen_cme, ycen_cme, $ pixrs, xloc, yloc, nloc, tscope) print, 'xm: ', xm aoi_image = TVRD () dotpos = STRPOS (cme_fits_file, '.fts') ; gifname = cme_fits_file ; STRPUT, gifname, '.gif', dotpos ; print, 'Output Image : ', gifname print, "AOI file: ", aoi_file TVLCT, red, green, blue, /GET WRITE_GIF, aoi_file, aoi_image, red, green, blue RETURN END