FUNCTION JMAP_CONSTRUCT2a, mvilist, angle, TIME0 = time0, TIME1 = time1, C2 = c2, $ C3 = c3, NOMOVIE = nomovie, PFILTER=pfilter, DFILTER=dfilter, RANGE = range, $ CHOP = dochop, FITSDIR=fitsdir, PROJ_TYPE=proj_type, SPHER=spher, width=width_pix, $ COORDS=coords, FRAME_SKIP=frame_skip, debug=debug, check_angs=check_angs ;+ ; $Id: jmap_construct2a.pro,v 1.26 2013/01/18 14:59:59 nathan Exp $ ; NAME: ; JMAP_CONSTRUCT2a ; ; PURPOSE: ; Perform actual jMap construction from a series of movies. ; ; CATEGORY: ; jMaps ; ; CALLING SEQUENCE: ; Result = JMAP_CONSTRUCT2a(mvilist, angle) ; ; INPUTS: ; mvilist: The list of movie files (.mvi). ; Angle: The position angle at which to make the map. ; ; OPTIONAL INPUTS: ; Time0: The starting time of the map. If not set, then the starting ; time of the first movie is used. ; Time1: The ending time of the map. If not set, then the ending time ; of the last movie is used. ; ; KEYWORD PARAMETERS: ; NOMOVIE: If set, the movie is not shown while the ; map is made. ; PFILTER: If set, applies a point filer to the map. ; DFILTER: If set, applies the standard deviation normalization ; filter to the map. ; CHOP: If set, chops the map. ; RANGE: desired range of the JMAP in degrees from sun center (start and end) ; ; OUTPUTS: ; Returns the jMap. ; ; RESTRICTIONS: ; In this version, if the starting and ending dates are not ; specified, then the movies should be in chronological order ; in the mvilist. ; ; PROCEDURE: ; Read the code. Note that the parameters for the optional ; image processing routines are constants found at the ; beginning of the program. ; ; EXAMPLE: ; Map = JMAP_CONSTRUCT2a(mvilist, angle, /C3) ; ; MODIFICATION HISTORY: ; Written by: Jeff Walters, Aug 1997 ; ; $Log: jmap_construct2a.pro,v $ ; Revision 1.26 2013/01/18 14:59:59 nathan ; move ishi sec_pix * 3600 ; ; Revision 1.25 2012/12/19 11:44:05 mcnutt ; replaced variable list with mvilist ; ; Revision 1.24 2012/11/14 18:13:35 nathan ; fix some debug stuff ; ; Revision 1.23 2012/06/20 18:01:59 nathan ; Do correction for incorrect xcen/ycen in early 2012 MVI files ; ; Revision 1.22 2012/05/22 16:53:58 mcnutt ; fix filename strput to work with mass _2m file names ; ; Revision 1.21 2011/03/23 19:14:01 nathan ; fix rectified-check logic ; ; Revision 1.20 2011/02/14 23:29:27 nathan ; Use same method(s) of determining attitude of MVI and instrument FOV for ; LASCO and SECCHI; verify ROLL of MVI file with user ; ; Revision 1.19 2011/02/08 19:50:17 nathan ; fix LASCO file search ; ; Revision 1.18 2010/07/14 19:35:33 sheeley ; nr - modify debug messages ; ; Revision 1.17 2009/11/13 21:03:48 nathan ; workaround for TB filenames in MVI header with extra char ; ; Revision 1.16 2009/07/27 15:57:54 cooper ; fixed wcs hdr for bad movies ; ; Revision 1.15 2009/07/15 20:58:45 cooper ; added cor2_point and hi_fix_pointer for better combined jmaps ; ; Revision 1.14 2009/07/15 16:51:32 cooper ; added check_angs keyword ; ; Revision 1.13 2009/04/03 18:02:36 nathan ; added debug messages ; ; Revision 1.12 2009/02/19 19:14:45 nathan ; check for subfields and verify ; ; Revision 1.11 2008/09/18 15:57:10 sheeley ; commit change from 9/12 (below) ; ; ; 080912, nbr - adjusted rotangle if FITS file not found ; ; Revision 1.10 2008/08/21 18:02:12 nathan ; implemented matchfitshdr2mvi for SECCHI only to properly account for roll and subfield of mvi ; ; Revision 1.9 2008/06/26 19:15:00 nathan ; added debug features ; ; Revision 1.8 2008/05/12 21:54:56 sheeley ; nr - revert back to r1.5 plus bug fix ; ; Revision 1.5 2008/03/31 15:56:47 nathan ; changed SOHO sun distance from .9 to .99 * 1au ; ; Revision 1.4 2008/02/06 18:57:51 nathan ; added DSUN_OBS keyword to plane_project.pro ; ; Revision 1.3 2007/11/08 18:22:17 nathan ; fix error from 11/7/07 edit ; ; Aug 2001 Jeff Walters. Current version. ; 2003.09.23, NRich - Add check for file_hdr.rectified ; 07.10.26, NRich - Added mvi version to sccmvihdr2struct calls ; ;- ;***************Notes****************** ;** If you want to use the start and end times of ;** the provided movies, then either do not set the ;** TIME0 and TIME1 keywords or set them to the value -1 ;** ;** Added keywords to control the use of point filters and ;** standard deviation filters. ;** Constants ;** PFILTER1_BOXSIZE = 3 PFILTER1_TOL = 7 PFILTER1_ITER = 2 PFILTER2_BOXSIZE = 3 PFILTER2_TOL = 7 PFILTER2_ITER = 2 DFACTOR = 1 ;always should be 1 CHOP_FACTOR = 5 IF(datatype(proj_type) EQ 'UND') THEN proj_type = 0 ; help,proj_type IF(~keyword_set(FRAME_SKIP)) THEN frame_skip = 1 spher = proj_type EQ 1 tai = long(0) n_images = 0 nfiles = n_elements(mvilist) ; hdr = def_mvi_strct(/secchi) ; this is nbr's change hdr = def_mvi_strct() ;nrs 10/25/07 ; update for new rev of dev_mvi_struct - nbr, 10/16/07 hdrs = REPLICATE(hdr,1) original_angle = angle n_ang = n_elements(angle) rotangle = replicate(-500.0, nfiles) ; rotangle is the CROTA value of the Movie (usually 0) sec_pix = fltarr(nfiles) rad_per_pix = fltarr(nfiles) sun_center = fltarr(nfiles, 2) proj = bytarr(nfiles) hi = bytarr(nfiles) xoffset = intarr(nfiles) size = intarr(nfiles, 2) frame_limit = intarr(nfiles, 2) Dsun = 0.99*1.5e11 ; SOHO distance to sun in meters width = 0 cols = 0 ;** Determine whether or not times have been specified by the user. user_times = 0 ;** Boolean- True means the user has specified ;** start and end times IF (KEYWORD_SET(time0) AND KEYWORD_SET(time1)) THEN BEGIN IF ( (time0 gt 0) AND (time1 gt 0) ) THEN user_times = 1 ENDIF ELSE BEGIN time0 = 1e+100 time1 = -1e+100 ENDELSE IF keyword_set(DEBUG) THEN BEGIN help, angle, TIME0, TIME1, C2, c3, NOMOVIE, PFILTER, DFILTER, RANGE, $ dochop, fitsdir, proj_type, spher, width_pix, coords, frame_skip help,mvilist print,mvilist ENDIF ;** Open file headers to determine image sizes. ;** and start/end times of movie FOR i=0, nfiles - 1 DO BEGIN OPENR,lu,mvilist[i],/GET_LUN SCCREAD_MVI, lu, file_hdr, ihdrs, imgs, swapflag rectified=file_hdr.rectified ;first_hdr = sccmvihdr2struct(ihdrs(0)) ;last_hdr = sccmvihdr2struct(ihdrs(file_hdr.nf-1)) first_hdr = sccmvihdr2struct(ihdrs(0),file_hdr.ver) mviroll=first_hdr.roll ; nbr, 10/26/07 - added mvi version last_hdr = sccmvihdr2struct(ihdrs(file_hdr.nf-1),file_hdr.ver) ;** Find start and end times if they are not already defined. frame_limit[i,0] = 0 frame_limit[i,1] = file_hdr.nf-1 frame_limit[i,1] = frame_limit[i,1] - (frame_limit[i,1] MOD frame_skip) curr_nim = (frame_limit[i,1] - frame_limit[i,0]) / frame_skip + 1 IF(NOT(user_times)) THEN BEGIN first_hdr = sccmvihdr2struct(ihdrs[frame_limit[i,0]],file_hdr.ver) ; nbr, 10/26/07 last_hdr = sccmvihdr2struct(ihdrs[frame_limit[i,1]],file_hdr.ver) ; nbr, 10/26/07 start_time = first_hdr.date_obs+' '+first_hdr.time_obs end_time = last_hdr.date_obs+' '+last_hdr.time_obs IF keyword_set(debug) then help,start_time,end_time t0 = LONG(ANYTIM2TAI(start_time)) t1 = LONG(ANYTIM2TAI(end_time)) IF(t0 LT time0) THEN time0 = t0 IF(t1 GT time1) THEN time1 = t1 FOR k=0,curr_nim-1 DO BEGIN j = frame_limit[i,0] + k*frame_skip img_hdr = sccmvihdr2struct(ihdrs[j],file_hdr.ver) ; nbr, 10/26/07 str = string(img_hdr.date_obs) + ' ' + string(img_hdr.time_obs) curr_tai = long(utc2tai(str2utc(str))) tai = [tai, curr_tai] ENDFOR ENDIF ELSE BEGIN FOR k=0,curr_nim-1 DO BEGIN ;j = frame_limit[i,0] + k*frame_skip;caryn commented this out 08/23/07 (you might need to change it back to the way it was at some point if this action caused error) j = k*frame_skip img_hdr = sccmvihdr2struct(ihdrs[j],file_hdr.ver) ; nbr, 10/26/07 str = STRING(img_hdr.date_obs) + ' ' + STRING(img_hdr.time_obs) IF keyword_set(debug) then help,str curr_tai = LONG(UTC2TAI(STR2UTC(str))) IF(curr_tai LT time0) THEN frame_limit[i,0] = j+frame_skip $ ELSE IF(curr_tai GT time1) THEN BEGIN frame_limit[i,1] = j-frame_skip BREAK ENDIF ELSE tai = [tai, curr_tai] ENDFOR IF(frame_limit[i,0] GE file_hdr.nf OR frame_limit[i,1] LT 0) THEN CONTINUE ; movie has no frames in chosen interval first_hdr = sccmvihdr2struct(ihdrs[frame_limit[i,0]],file_hdr.ver) ; nbr, 10/26/07 last_hdr = sccmvihdr2struct(ihdrs[frame_limit[i,1]],file_hdr.ver) ; nbr, 10/26/07 ENDELSE curr_nim = (frame_limit[i,1] - frame_limit[i,0]) / frame_skip + 1 found = 0 print,'Checking ',first_hdr.filename fname='' tel = STRMID(first_hdr.filename, 18,2) IF (tel EQ 'h1' OR tel EQ 'h2') THEN hi[i] = 1 ishi=hi[i] iscor2=(tel EQ 'c2') IF(STRLEN(first_hdr.filename) EQ 12 AND proj_type GT 0) THEN BEGIN ; LASCO tel = first_hdr.detector ;IF(tel EQ 'C2' OR tel EQ 'C3') THEN lasco[i] = 1 IF keyword_set(FITSDIR) THEN fname = concat_dir(fitsdir,first_hdr.filename) $ ELSE BEGIN ; first check LZ fname=lz_from_ql(first_hdr) ; then check QL IF fname EQ 'none' THEN fname=lz_from_ql(first_hdr,/ql) ENDELSE IF file_exist(fname) THEN BEGIN found = 1 im = lasco_readfits(fname, hdr0,/no_img) scchdr=convert2secchi(hdr0) ENDIF ELSE message,'ERROR: '+fname+' not found.' ENDIF ELSE $ IF(STRLEN(first_hdr.filename) EQ 25) THEN BEGIN ; SECCHI ;IF(NOT(keyword_set(FITSDIR))) THEN fname = scc_findfile(first_hdr.filename) $ ; some HI2 movies have weird values for ccd_pos; do double-check ; based on R1COL value osc='A' IF strmid(first_hdr.filename,strpos(first_hdr.filename,'.f')-1,1) NE 'A' THEN osc='B' ccdpos=file_hdr.ccd_pos IF total(ccdpos) NE 0 THEN BEGIN print,'R1ROW, R1COL, R2ROW, R2COL:',ccdpos IF (ishi and osc EQ 'A' and ccdpos[1] NE 51) OR $ (ishi and osc EQ 'B' and ccdpos[1] NE 79) THEN BEGIN inp='' print,'' print,'If this is a subfield, enter y' read, 'else return if this is a FFV movie: ',inp IF inp EQ '' THEN file_hdr.ccd_pos[*]=0 ENDIF ENDIF ; end ccdpos check IF(NOT(keyword_set(FITSDIR))) THEN BEGIN filen=first_hdr.filename fname = sccfindfits(filen) IF fname EQ '' THEN BEGIN ; nbr, 11/7/07 - workaround for filename from secchi_prep calibrated image strput,filen,'??',16 fname = sccfindfits(filen) help,fname ENDIF IF fname EQ '' THEN BEGIN ; workaround for TB filenames break_file,filen,di,pa,fi,su last4=rstrmid(fi,4) filen=strmid(fi,0,16)+'?'+last4+'.fts' fname=sccfindfits(filen) help,fname ENDIF ENDIF ELSE fname = fitsdir + '/' + first_hdr.filename IF(fname NE '') THEN BEGIN found = 1 im = sccreadfits(fname, scchdr, /NODATA) ENDIF ENDIF IF (found) THEN BEGIN ; SECCHI or LASCO crota=scchdr.crota IF rectified GT 180 THEN rectified=rectified-360. IF ABS(rectified + crota) GT 1 and mviroll EQ 0 THEN BEGIN window, 0, xsize = 512, ysize = 512 tv,imgs[0] print,'' help,crota,rectified,mviroll print,'Cannot tell from MVI header if North is up in MVI.' inp=0 print,' Please enter the correct value for CROTA of the MVI frames ' read,' (Enter zero if you believe Solar North is up): ',inp mviroll = float(inp) ;clockwise correction angle file_hdr.rectified=round(mviroll-crota) first_hdr.roll=mviroll ENDIF help,mviroll rotangle[i]=-mviroll IF(ishi) THEN BEGIN hi_fix_pointing,scchdr ; hi_fix_pointing changes CRVAL, not CRPIX. first_hdr.xcen=0 first_hdr.ycen=0 ; Do correction for incorrect xcen/ycen in early 2012 MVI files. ; CRPIX set to image center in matchfitshdr2mvi.pro. ENDIF IF(iscor2) THEN cor2_point,scchdr,/nojitter,/rollzero ; Correct roll/subfield is in mvi file, but must also have header matchfitshdr2mvi,scchdr,file_hdr,first_hdr ; this accounts for changes in cdelt, roll, subfield dsun=scchdr.dsun_obs sec_pix[i] = scchdr.cdelt1 ;spanx = file_hdr.ccd_pos[2] - file_hdr.ccd_pos[0] ;IF(spanx GT 0) THEN sec_pix[i] = sec_pix[i] * spanx / float(scchdr.r2row-scchdr.r1row) IF(proj_type GT 0) and (rectified LT 10) THEN BEGIN ; make sure image has not been rotated 180 deg wcs = FITSHEAD2WCS(scchdr) IF(scchdr.detector EQ 'HI1') THEN BEGIN ;correction in wcshdr for movies made in late 07, early 08 IF(ABS(wcs.cdelt[0]) LT .001) THEN BEGIN ;tjc 2009-07-27 wcs.cdelt=[0.0399,0.0399] wcs.crpix=[256.5,256.5] ENDIF ENDIF IF(scchdr.detector EQ 'HI2') THEN BEGIN IF(ABS(wcs.cdelt[0]) LT .001) THEN BEGIN wcs.cdelt=[.144,.144] wcs.crpix=[256.5,256.5] ENDIF ENDIF ; NOTE: This is all done in matchfitshdr2mvi above: ;x1 = 0 & x2 = 0 & y1 = 0 & y2 = 0 ;IF(file_hdr.ccd_pos[1] NE 0) THEN BEGIN ; spanx = float(scchdr.r2row - scchdr.r1row) ; spany = float(scchdr.r2col - scchdr.r1col) ; x1 = scchdr.naxis1 * float(file_hdr.ccd_pos[0]-scchdr.r1row) / spanx ; x2 = scchdr.naxis1 * float(file_hdr.ccd_pos[2]-scchdr.r1row) / spanx ; y1 = scchdr.naxis2 * float(file_hdr.ccd_pos[1]-scchdr.r1col) / spany ; y2 = scchdr.naxis2 * float(file_hdr.ccd_pos[3]-scchdr.r1col) / spany ; IF keyword_set(DEBUG) THEN help,x1,x2,y1,y2 ;ENDIF proj[i] = 1 coord = WCS_GET_COORD(wcs) ; array of coordinates for all pixels in original image ;sz = size(coord) ;IF(file_hdr.ccd_pos[1] NE 0) THEN BEGIN ; movie is from cropped images ; IF(x1 GE 0 AND x2 GE 0 AND y1 GE 0 AND y2 GE 0 AND $ ; x1 LT sz[1] AND x2 LT sz[1] AND y1 LT sz[2] AND y2 LT sz[2]) THEN $ ; coord = (temporary(coord))[*, x1:x2-1, y1:y2-1] ;ENDIF IF(~hi[i]) THEN coord = temporary(coord) / 3600. ; others are in arcseconds ENDIF ELSE BEGIN print,'*** NOT using spherical projection for ',scchdr.detector,' ***' ;reg = (file_hdr.ccd_pos[1] NE 0) ? [x2-x1, y2-y1] : wcs.NAXIS ;sun_center[i,*] = (WCS_GET_PIXEL(wcs,[0,0]) - [x1,y1]) * float([file_hdr.nx, file_hdr.ny]) / reg ;sec_pix[i] = scchdr.CDELT1 * (float(scchdr.NAXIS1) / file_hdr.nx) ;sun_center[i,*] = WCS_GET_PIXEL(wcs,[0,0]) sun_center[i,*]=[file_hdr.sunxcen,file_hdr.sunycen] ENDELSE ENDIF ; NOW figure out roll of MVI frames IF keyword_set(DEBUG) THEN BEGIN print,'sec_pix=',sec_pix print,'sun_center=',sun_center help,crota,mviroll,rectified wait,2 device,get_screen_size=scrsiz ENDIF IF ~found THEN BEGIN sec_pix[i] = file_hdr.sec_pix sun_center[i,0] = file_hdr.sunxcen sun_center[i,1] = file_hdr.sunycen ;nbr, 9/12/08 rotangle[i]=-mviroll print, 'jmap_construct2: Could not find source FITS file--using mvi header values' ENDIF IF (ishi) THEN sec_pix[i] = sec_pix[i] * 3600. ;HI scales are in degrees hdrs = [hdrs, first_hdr] IF(size[i,0] EQ 0) THEN size[i,*] = [file_hdr.nx, file_hdr.ny] IF(NOT(KEYWORD_SET(nomovie))) THEN $ window, 0, xsize = 512, ysize = 512, title = 'jMaps: Movie' IF(proj[i]) THEN BEGIN ; project in blocks to avoid memory allocation failure proj_block = floor(50.*512.*512. / (float(file_hdr.nx > 512) * float(file_hdr.ny > 512))) proj_count = 0 proj_nim = curr_nim < proj_block new_imgs = bytarr(proj_nim, file_hdr.nx, file_hdr.ny) FOR j=0,proj_nim-1 DO new_imgs[j,*,*] = imgs[j*frame_skip + frame_limit[i,0]] print, '% JMAP_CONSTRUCT2A: projecting images', 1, ' to', proj_nim, ' for '+mvilist[i] IF(datatype(proj_map) NE 'UND') THEN undefine, proj_map PLANE_PROJECT, new_imgs, coord, proj_imgs, LIMIT=range[1], SPHER=spher, SCEN=scen, RANGE=rad, $ PIXEL_MAP=proj_map, /INTERP, NEW_COORD=proj_coord, DSUN_OBS=dsun, debug=debug s = size(proj_imgs) size[i,*] = s[2:3] rad_per_pix[i] = (rad[1]-rad[0]) / s[2] range[1] = max(abs(rad)) < range[1] sun_center[i,*] = scen rotangle[i] = 0. ENDIF ELSE print,'% JMAP_CONSTRUCT2A: **NOT** projecting images for '+mvilist[i] IF(rotangle[i] LT (-400.)) THEN BEGIN ; this means rotangle still has init value. sun_center[i,*] = [file_hdr.sunxcen, file_hdr.sunycen] IF (file_hdr.rectified GT 0) THEN rotangle[i]=0 ELSE BEGIN ; NRich, 03.09.23 pnt = GET_SC_POINT(str) rotangle[i] = -1*pnt.sc_roll ; degrees message,'Using '+trim(-1*pnt.sc_roll)+' from GET_SC_POINT for roll angle.',/info ENDELSE IF(sun_center[i,0] EQ 0) THEN BEGIN MESSAGE, 'Incomplete file header. Calculating sun center.', /INFORMATIONAL fitshdr = STRUCT2FITSHDR(img_hdr) sun = GET_SUN_CENTER(fitshdr, FULL = file_hdr.nx) file_hdr.sunxcen = sun.xcen file_hdr.sunycen = sun.ycen MESSAGE, 'Sun Center = ['+ STRTRIM(file_hdr.sunxcen) + ', '+STRTRIM(file_hdr.sunycen)+']', /INFORMATIONAL ENDIF ENDIF IF(~proj[i]) THEN rad_per_pix[i] = sec_pix[i] / 3600. curr_cols = floor((range[1] - range[0]) / rad_per_pix[i]) xoffset[i] = (range[0] / rad_per_pix[i]) curr_width=round(width_pix/rad_per_pix[i]) IF ((floor(curr_width/2)-(curr_width/2.)) EQ 0) THEN curr_width=curr_width+1 IF(curr_cols GT cols) THEN cols = curr_cols IF(curr_width GT width) THEN width = curr_width width_pix=width*rad_per_pix[i] widthp = width/2 + 1 angle = original_angle + 90 + rotangle[i] IF keyword_set(DEBUG) THEN help, angle,width,cols,curr_cols, curr_width, curr_nim, n_ang curr_strips = bytarr(curr_cols, curr_width, curr_nim, n_ang) s=size[i,*] hyp = ceil(sqrt(LONG(s[0])^2 + LONG(s[1])^2)) rot_image = BYTARR(hyp,hyp) rot_cen = fltarr(2, n_ang) rot_matrix = fltarr(2, 2, n_ang) rot_pos = fltarr(5,2) scen = reform(sun_center[i,*]) - [s[0]/2, s[1]/2] + [hyp/2, hyp/2] FOR k=0,n_ang-1 DO BEGIN ang = angle[k] * !PI/180. rot_matrix[*,*,k] = [[cos(ang), -sin(ang)], [sin(ang), cos(ang)]] ; rotates CW by angle rot_cen[*,k] = (rot_matrix[*,*,k] # (reform(sun_center[i,*]) - [s[0]/2, s[1]/2])) + [hyp/2, hyp/2] ENDFOR IF(proj[i]) THEN BEGIN ; store coordinates along strip axis coord_image = fltarr(2,hyp,hyp) coords = fltarr(2,cols) FOR j=0,1 DO BEGIN ; coord_image[j, (hyp-s[0])/2:(hyp-s[0])/2 + s[0]-1, (hyp-s[1])/2:(hyp-s[1])/2 + s[1]-1] = proj_coord[j,*,*] ; coord_image[j,*,*] = rot(reform(coord_image[j,*,*]), angle[0]) ; coords[j,*] = EXTRAC(reform(coord_image[j,*,*]), xcen, ycen, curr_cols, 1) ENDFOR ENDIF ELSE message,'Not defining coords ',/info ;IF keyword_set(DEBUG) THEN window,xsize=(hyp