;+ ; fmc.pro IDL procedure ; ; PURPOSE ; Compute absolute/excess coronal mass. ; Read FITS images and define a region of interest. ; NOTE: Works for LASCO C1/C2/C3 images only. ; ; SYNTAX ; fmc, 'cme_fits_file' ; fmc, 'cme_fits_file', 'ref_fits_file' ; ; PROCEDURES USED ; fxpar Get FITS parameter associated with a keyword ; defroic (tailored version of "defroi") ; get_cal_name ; readfits Read FITS file ; write_gif Write GIF file ; kb_cyl Keyboard input to define curved wedge ; mouse_cyl Mouse cursor to define curved wedge ; kb_wedge Keyboard input to define straight wedge ; mouse_wedge Mouse cursor to define straight wedge ; mouse_poly Mouse cursor to define polygon ; file_poly File to define polygon ; xmass compute excess/absolute mass ; ; HISTORY ; Author: Andrew L. Stanger 2 June 1998 ; 1 Oct 1998: Add logic for C1, C2, C3, EIT telescopes. ; 2 Oct 1998: Use 'dte' & 'get_cal_name' for LASCO telescopes ONLY. ; 19 Oct 1998: Add scope_name argument to xmass function. ; Date order is reversed from FITS convention for DATE-OBS. ; 10 Nov 1998: Change 'defroi' name to 'defroic' to distinguish ; it from the RSI version. ; 28 Jan 1999: version 3: Add ability to store/fetch polygon vertex file. ; 29 Jan 1999: version 4: Add ability to define a wedge area. ; 8 Feb 1999: version 5: Add ability to define a cylindrical wedge area. ; 11 Feb 1999: Use FINDFILE to avoid overwriting *.pv, *.gif, *.fm files. ; 12 Feb 1999: Make cmeyear 4 digits, not 2 (i.e., Y2K compliance). ;- PRO fmc, cme_fits_file, ref_fits_file if (N_PARAMS() NE 2) THEN ref_fits_file = '' ;--- Initialize variables. i = 0L nloc = 0L regnum = 0 red = bytarr (256) green = bytarr (256) blue = bytarr (256) mx = 0 ; =0 cme image only, =1 cme & ref images tscope = 0 ; Telescope designation ; 1: C1 2: C2 3: C3 4: EIT 5: SXT ; 6: MK3 7: SMM C/P 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 hdu_cme = '' ; CME header telescop_cme = '' ; Telescope instrume_cme = '' ; Instrument detector_cme = '' ; Detector roll_cme = 0.0 ; S/C roll. 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 hdu_ref = '' ; REF header telescop_ref = '' ; Telescope instrume_ref = '' ; Instrument detector_ref = '' ; Detector ;--- Read CME FITS image. ; imgcme = rfits (cme_fits_file, header=hdu_cme, /scale) imgcme = readfits (cme_fits_file, hdu_cme) dspcme = imgcme imin = MIN (imgcme, max=imax) print, 'imgcme: imin, imax: ', imin, imax ;--- Retrieve parameters from FITS header. 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") cdelt1_cme = fxpar (hdu_cme, "CDELT1") crradius_cme = fxpar (hdu_cme, "CRRADIUS") pixrs = crradius_cme crota1_cme = fxpar (hdu_cme, "CROTA1") scroll_cme = fxpar (hdu_cme, "SCROLL") roll_cme = -crota1_cme + scroll_cme print, 'fm> crradius_cme: ', crradius_cme, ' 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") instrume_cme = fxpar (hdu_cme, "INSTRUME") detector_cme = fxpar (hdu_cme, "DETECTOR") cmedate = fxpar (hdu_cme, "DATE-OBS") cmetime = fxpar (hdu_cme, "TIME-OBS") print, 'cmedate: ', cmedate, ' STRLEN (cmedate): ', STRLEN (cmedate) print, 'cmetime: ', cmetime ;--- Extract day, month & year from 'DATE-OBS'. IF (tscope GT 3) THEN $ BEGIN cmeday = STRMID (cmedate, 0, 2) cmemonth = STRMID (cmedate, 3, 2) cmeyear = STRMID (cmedate, 6, STRLEN (cmedate) - 6) print, 'STRLEN (cmedate): ', STRLEN (cmedate) print, 'cmeyear:...', cmeyear, '...' IF (STRLEN (cmeyear) EQ 2) THEN $ BEGIN IF (cmeyear GE 50) THEN cmeyear = cmeyear + 1900 $ ELSE cmeyear = cmeyear + 2000 END print, 'STRLEN (cmeyear): ', STRLEN (cmeyear) IF (STRLEN (cmeyear) GT 4) THEN $ cmeyear = STRCOMPRESS (STRTRIM (STRMID (cmeyear,STRLEN (cmeyear)-4,4))) print, 'cmeyear:...', cmeyear, '...' yymmdd = cmeyear + cmemonth + cmeday print, 'fm> yymmdd: ', yymmdd END $ ;--- LASCO/C1/C2/C3: ELSE $ ; LASCO BEGIN dte = STR2UTC (cmedate) yymmdd = UTC2YYMMDD (dte) print, 'yymmdd: ', yymmdd, ' STRLEN (yymmdd): ', STRLEN (yymmdd) print, 'size (yymmdd): ', size (yymmdd) IF (STRLEN (yymmdd) EQ 6) THEN $ BEGIN cmeyear = STRMID (yymmdd, 0, 2) cmemonth = STRMID (yymmdd, 2, 2) cmeday = STRMID (yymmdd, 4, 2) IF (cmeyear GE 50) THEN cmeyear = cmeyear + 1900 $ ELSE cmeyear = cmeyear + 2000 END $ ELSE $ BEGIN cmeyear = STRMID (yymmdd, 0, 4) cmemonth = STRMID (yymmdd, 4, 2) cmeday = STRMID (yymmdd, 6, 2) END cmeyear = STRMID (cmeyear, STRLEN (cmeyear) - 4, 4) print, 'STRLEN (cmeyear): ', STRLEN (cmeyear) print, 'cmeyear:...', STRTRIM (STRING (cmeyear)), '...' print, 'STRLEN (cmeyear): ', STRLEN (cmeyear) END ;--- Extract time from TIME-OBS. cmehour = STRMID (cmetime, 0, 2) cmeminute = STRMID (cmetime, 3, 2) cmesecond = STRMID (cmetime, 6, 2) print, 'cmeyear: ', cmeyear print, 'cmemonth: ', cmemonth print, 'cmeday: ', cmeday PRINT, 'cme year, month, day: ', cmeyear, cmemonth, cmeday PRINT, 'cme hour, minute, second: ', cmehour, cmeminute, cmesecond ;--- Perform telescope-specific operations as needed. ;--- Solar Maximum Mission Coronagraph/Polarimeter: IF (STRPOS (telescop_cme, "SMM C/P") NE -1) THEN $ ; SMM C/P image. BEGIN tscope = 7 dspcme = imgcme ; END ;--- Mauna Loa Mark-III K-coronameter: IF (STRPOS (telescop_cme, "MK3") NE -1) THEN $ ; MK3 image. BEGIN ;{ tscope = 6 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 ;} IF (STRPOS (telescop_cme, "SOHO") NE -1) THEN $ ; SOHO BEGIN ;{ IF (STRPOS (instrume_cme, "LASCO") NE -1) THEN $ ; LASCO BEGIN ;{ ;--- SOHO C1 Telescope: IF (STRPOS (detector_cme, "C1") NE -1) THEN $ ; C1 BEGIN ;{ tscope = 1 END $ ;} ;--- SOHO C2 Telescope: ELSE IF (STRPOS (detector_cme, "C2") NE -1) THEN $ ; C2 BEGIN ;{ tscope = 2 sd = GETENV ('NRL_LIB') + '/lasco/data/calib/' ; dte = STR2UTC (cmedate) ; yymmdd = UTC2YYMMDD (dte) vig_fn = get_cal_name (sd + 'C2_cl*vig*.dat', yymmdd) IF (vig_fn NE '') THEN $ BEGIN ;{ vig = readfits (vig_fn, vig_hdu) dspcme = imgcme / vig ENDIF $ ;} ELSE BEGIN ;{ print, 'ERROR: fmc - No ' + sd + 'C2_cl*vig*.dat file' dspcme = imgcme * 0.0 RETURN END ;} END $ ;} ;--- SOHO C3 coronagraph: ELSE IF (STRPOS (detector_cme, "C3") NE -1) THEN $ ; C3 BEGIN ;{ tscope = 3 sd = GETENV ('NRL_LIB') + '/lasco/data/calib/' vig_fn = get_cal_name (sd + 'C3_cl*vig*.dat', yymmdd) IF (vig_fn NE '') THEN $ BEGIN ;{ vig = FLTARR (1024,1024) OPENR, lucal, vig_fn, /GET_LUN READU, lucal, vig ; IF (szlog (0) NE 0) THEN PRINTF, lulog, 'Used ' + vig_fn PRINT, 'Used ' + vig_fn CLOSE, lucal FREE_LUN, lucal dspcme = imgcme / vig ENDIF $;} ELSE BEGIN ;{ print, 'ERROR: c3_calibrate - No ' + sd + 'C3_cl*vig*.dat file' dspcme = imgcme * 0.0 RETURN END ;} END ;} END $ ;} ;--- SOHO Extreme Ultra-Violet Imaging Telescope: ELSE IF (STRPOS (instrume_cme, "EIT") NE -1) THEN $ ; EIT BEGIN ;{ tscope = 4 END ;} END ;} ;--- Yohkoh Soft X-Ray Telescope: IF (STRPOS (telescop_cme, "SXT") NE -1) THEN $ ; SXT BEGIN tscope = 5 END ;--- Use NRL_LIB version of pb0r: date_time = cmedate + ' ' + cmetime solar_info = pb0r (date_time) ; Get solar P, B0, Rsun. Pangle = solar_info [0] ; P-angle B0 = solar_info [1] ; B0 angle Rsunarc = solar_info [2] * 60.0 ; Solar radius (arcsec). pixrs = Rsunarc / cdelt1_cme ; Pixels / Rsun. print, 'fmc> pixrs: ', pixrs 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) imgref = readfits (ref_fits_file, hdu_ref) 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") instrume_ref = fxpar (hdu_ref, "INSTRUME") detector_ref = fxpar (hdu_ref, "DETECTOR") refdate = fxpar (hdu_ref, "DATE-OBS") reftime = fxpar (hdu_ref, "TIME-OBS") IF (tscope GT 3) THEN $ BEGIN refday = STRMID (refdate, 0, 2) refmonth = STRMID (refdate, 3, 2) refyear = STRMID (refdate, 6, STRLEN (refdate) - 6) IF (STRLEN (refyear) EQ 2) THEN $ BEGIN IF (refyear GE 50) THEN refyear = refyear + 1900 $ ELSE refyear = refyear + 2000 END yymmdd = refyear + refmonth + refday END $ ;--- LASCO/C1/C2/C3: ELSE $ BEGIN dte = STR2UTC (refdate) yymmdd = UTC2YYMMDD (dte) IF (STRLEN (yymmdd) EQ 6) THEN $ BEGIN refyear = STRMID (yymmdd, 0, 2) refmonth = STRMID (yymmdd, 2, 2) refday = STRMID (yymmdd, 4, 2) IF (refyear GE 50) THEN refyear = refyear + 1900 $ ELSE refyear = refyear + 2000 END $ ELSE $ BEGIN refyear = STRMID (yymmdd, 0, 4) refmonth = STRMID (yymmdd, 4, 2) refday = STRMID (yymmdd, 6, 2) END END refhour = STRMID (reftime, 0, 2) refminute = STRMID (reftime, 3, 2) refsecond = STRMID (reftime, 6, 2) base_name = STRTRIM (STRING (cmeyear)) $ + STRING (cmemonth) + STRING (cmeday) $ + '.' $ + STRING (cmehour) + STRING (cmeminute) + '-' $ + STRING (refhour) + STRING (refminute) 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 = 7 dspref = imgref END if (STRPOS (telescop_ref, "MK3") NE -1) THEN $ ; MK3 image. BEGIN tscope = 6 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 IF (STRPOS (telescop_ref, "SOHO") NE -1) THEN $ ; SOHO BEGIN ;{ IF (STRPOS (instrume_ref, "LASCO") NE -1) THEN $ ; LASCO BEGIN ;{ IF (STRPOS (detector_ref, "C1") NE -1) THEN $ ; C1 BEGIN ;{ tscope = 1 END $ ;} ELSE IF (STRPOS (detector_ref, "C2") NE -1) THEN $ ; C2 BEGIN ;{ tscope = 2 sd = GETENV ('NRL_LIB') + '/lasco/data/calib/' ; dte = STR2UTC (cmedate) ; yymmdd = UTC2YYMMDD (dte) vig_fn = get_cal_name (sd + 'C2_cl*vig*.dat', yymmdd) IF (vig_fn NE '') THEN $ BEGIN ;{ vig = readfits (vig_fn, ref_hdu) PRINT, 'Used ' + vig_fn dspref = imgref / vig ENDIF $ ;} ELSE BEGIN ;{ PRINT,'ERROR: fmc - No ' + sd + 'C2_cl*vig*.dat file' RETURN END ;} END $ ;} ELSE IF (STRPOS (detector_ref, "C3") NE -1) THEN $ ; C3 BEGIN ;{ tscope = 3 sd = GETENV ('NRL_LIB') + '/lasco/data/calib/' ; dte = STR2UTC (cmedate) ; yymmdd = UTC2YYMMDD (dte) ; print, 'fm> dte: ', dte, ' yymmdd: ', yymmdd vig_fn = get_cal_name (sd + 'C3_cl*vig*.dat', yymmdd) IF (vig_fn NE '') THEN $ BEGIN ;{ vig = FLTARR (1024,1024) OPENR, lucal, vig_fn, /GET_LUN READU, lucal, vig PRINT, 'Used ' + vig_fn CLOSE, lucal FREE_LUN, lucal dspref = imgref / vig ENDIF $;} ELSE BEGIN ;{ print, 'ERROR: fmc - No ' + sd + 'C3_cl*vig*.dat file' dspref = imgref * 0.0 RETURN END ;} END ;} END $ ;} ELSE IF (STRPOS (instrume_ref, "EIT") NE -1) THEN $ ; EIT BEGIN ;{ tscope = 4 END ;} END ;} IF (STRPOS (telescop_ref, "SXT") NE -1) THEN $ ; SXT BEGIN tscope = 5 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 7) THEN imax = imax * 0.5 IF (tscope EQ 6) THEN imax = imax * 0.2 IF (ABS (imin) LT abs (imax)) THEN imax = - imin IF (ABS (imin) GT abs (imax)) THEN imin = - imax print, 'dif min/max: ', imin, imax IF (tscope EQ 1) THEN $ ; C1 BEGIN imin = -5.0e-9 imax = 5.0e-9 END IF (tscope EQ 2) THEN $ ; C2 BEGIN imin = -5.0e-9 imax = 5.0e-9 END IF (tscope EQ 3) THEN $ ; C3 BEGIN imin = -5.0e-9 imax = 5.0e-9 END 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 base_name = STRTRIM (STRING (cmeyear)) $ + STRING (cmemonth) + STRING (cmeday) + '.' $ + STRING (cmehour) + STRING (cmeminute) END ;} CASE tscope OF 1: scope_name = 'soho_c1' 2: scope_name = 'soho_c2' 3: scope_name = 'soho_c3' 4: scope_name = 'soho_eit' 5: scope_name = 'yohkoh_sxt' 6: scope_name = 'mlso_mk3' 7: scope_name = 'smm_cp' ELSE: scope_name = 'Unknown' ENDCASE ;--- Action Loop loop: print, 'Request Action' IF (mx EQ 0) THEN TV, dspcmeb $ ; Display CME image. ELSE TV, dspdif ; Display difference image. regnum = regnum + 1 ; Increment region number. print, 'regnum: ', regnum print, 'base_name: ', base_name ;--- Determine valid region number & define output file names. outfiles: IF (regnum LE 9) THEN $ BEGIN pv_file = scope_name + '.' + base_name + '.r0' + $ STRTRIM (STRING (regnum, FORMAT='(I1)')) + '.pv' aoi_file = scope_name + '.' + base_name + '.r0' + $ STRTRIM (STRING (regnum, FORMAT='(I1)')) + '.gif' log_file = scope_name + '.' + base_name + '.r0' + $ STRTRIM (STRING (regnum, FORMAT='(I1)')) + '.fm' END $ ELSE $ BEGIN pv_file = scope_name + '.' + base_name + '.r' + $ STRTRIM (STRING (regnum, FORMAT='(I2)')) + '.pv' aoi_file = scope_name + '.' + base_name + '.r' + $ STRTRIM (STRING (regnum, FORMAT='(I2)')) + '.gif' log_file = scope_name + '.' + base_name + '.r' + $ STRTRIM (STRING (regnum, FORMAT='(I2)')) + '.fm' END ;--- Check for pre-existing mass log file. ; If it exists, increment the region number and try again. fname = FINDFILE (log_file, COUNT=nfiles) if (nfiles > 0) THEN $ BEGIN regnum = regnum + 1 GOTO, outfiles END print, 'fm> pv_file: ', pv_file print, 'fm> aoi_file: ', aoi_file print, 'fm> log_file: ', log_file ;--- Request Action Command. print, 'Refresh Display [d]' print, 'Keyboard input to define curved wedge [kc]' print, 'Mouse cursor to define curved wedge [mc]' print, 'Keyboard input to define wedge [kw]' print, 'Mouse cursor to define wedge [mw]' print, 'Mouse cursor to define polygon [p]' print, 'Read disk file to define polygon [f]' print, 'Quit program [q]' action = '' read, 'action: ', action CASE action OF 'd': BEGIN ; IF (mx EQ 0) THEN TV, dspcmeb $ ; Display CME image. ; ELSE TV, dspdif ; Display difference image. regnum = regnum - 1 GOTO, loop END 'kc': BEGIN kb_cyl5, xdim_cme, ydim_cme, xcen_cme, ycen_cme, pixrs, roll_cme, $ regnum, log_file, pv_file, xv, yv, nver, polypix END 'mc': BEGIN mouse_cyl5, xdim_cme, ydim_cme, xcen_cme, ycen_cme, pixrs, roll_cme, $ regnum, log_file, pv_file, xv, yv, nver, polypix END 'kw': BEGIN kb_wedge5, xdim_cme, ydim_cme, xcen_cme, ycen_cme, pixrs, roll_cme, $ regnum, log_file, pv_file, xv, yv, nver, polypix END 'mw': BEGIN mouse_wedge5, xdim_cme, ydim_cme, xcen_cme, ycen_cme, pixrs, roll_cme, $ regnum, log_file, pv_file, xv, yv, nver, polypix END 'p': BEGIN mouse_poly5, xdim_cme, ydim_cme, $ regnum, log_file, pv_file, xv, yv, nver, polypix END 'f': BEGIN file_poly5, xdim_cme, ydim_cme, $ regnum, log_file, pv_file, xv, yv, nver, polypix END 'q': RETURN ELSE: RETURN ENDCASE ;--- Determine coordinates of all pixel within polygon area. 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) nver = size_xv [1] ; print, 'size [xv]: ', size (xv) ; print, 'size [yv]: ', size (yv) ; FOR i = 0L, nver - 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 = xmass5 (log_file, scope_name, 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 GOTO, loop RETURN END