;+ ;NAME: ivm_mdi_align ; ;PURPOSE: fix the IVM pointing by using MDI mag, Ic images ; ;CATEGORY: data processing ; ;CALLING SEQUENCE: ivm_align_with_mdi, mag, oblk, newlat, newcmd, newxy ; ;INPUTS: mag: UH-type magnetic field structure, e.g. from mag=ivm_bfits(file), ; oblk: IVM header *structure* ; mask: Puka/limb mask - use to avoid problem pixels when computing stats ; ;OPTIONAL KEYWORD PARAMETERS: ; err - String containing error flag(s) ; blong - array containing blong - ELW-20100312 ; icont - Continuum intensity ; noflag - Do not log (with savegifs) if there was an error; assumably ; the sign is wrong and will be called again ; ref_offset - Reference pixel offset due to crop/box shift ; saveDir - Directory to save log file to (if log is set) ; suffix - Suffix to append to name in log ; log - If set, log the results = append to file in directory saveDir ; verbose - Verbosity [0-3] - if 3, will log the results as with /log ; ;OUTPUTS: ; NEWLAT: the adjusted latitude, in radians ; NEWCMD: the adjusted longitude, in radians ; ;COMMON BLOCKS: ;SIDE EFFECTS: ; ;RESTRICTIONS: MDI data files need to be available, location is hard-wired. ; MDI data level 1.8 required ; ;PROCEDURE: This program takes an ivm oblk and blong_qnd and uses the closest ; available mdi magnetogram file to correct the IVM pointing information ; ;METHOD: ; [NOT:] A piece of the IVM image is first selected. ; lookForSunspots draws a box around the likely region of the spot ; (where >.80 % abs(max(image)) is met). ; If FULL keyword is set, the entire IVM image is used, though masked ; areas (pukas, problem edges) are set to zero. Otherwise, the image ; is cropped to a smaller one. ; A larger, mdi magnetogram is used for alignment with the same central ; point as the ivm image. The function correl optimize is used to ; find the ideal x/y shift to align the images, factors out the shift ; due to the different sized images. This pixel shift is is translated ; back into the proper lat/lon offset. ; If mdi magnetogram is not available, alignment will be tried with ; closest mdi continuum file for that date. Cont. files typically work, ; but not when field strength is low, so magnetograms are preferred. ; ;MODIFICATION HISTORY: ; IVM_align_with_MDI.pro: ; TDunn, 2009, adopting some of TRM's automatic-alignment codes, ; adding automatic MDI-data-finding. ; KDL Sep 2009: standardizing various, and ensuring that all updating is ; done as advertized. ; KDL Oct 2009: removing interactive bits, and correcting terrestrial/L-1 ; offsets, following more closely TRM's ivm_mdi_align.pro but w/ this ; faster optimization. ; ivm_mdi_align.pro: ; ELW 12 Mar 2010: Renamed as this is for use with rivm, pre-inversion. ; Got rid of gotos; added Blong array keyword. ; ELW May 2010: No more "mag" -- again, this is rivm-centric. ; ELW Jun 2010: Shift scales no longer hard-coded ; ELW Oct 2013: Handle the case where MDI files contain >50% NaN ;- PRO ivm_mdi_align, oblk, newlat, newcmd, newxy, blong, mask, icont=icont, $ verbose=verbose, err=err, ref_offset=ref_offset, $ saveDir=saveDir, suffix=suffix, noflag=noflag, log=log if (keyword_set(verbose)) then vb=verbose else vb=0 ; Location/path information for images, logs, ouput ; Default saveDir if it is not passed in if ~keyword_set(log) then log=0 else log=1 if ~keyword_set(saveDir) then saveDir='/data/cora/wagneric/RivmTest/' ;Originally was: ;saveDir='/data/cora/SOLMAG/ANALYSIS/IVM_REDUCTION/QUICKLOOK_SURVEY/POINTING/' imageDir=saveDir+'gifs/' ;directory for saving images logFileDir=saveDir+'LOGS/' ; directory for log files of pointing changes err='' use_blong=1 pntng_chng=1 if (keyword_set(mag)) then usemag=1 else usemag=0 date_obs = ANYTIM( oblk.date_obs, /STIME, /DATE_ONLY ) time_obs = ANYTIM( oblk.date_obs, /STIME, /TIME_ONLY ) ; Note: not rounding off seconds ivmdate0 = anytim(date_obs+' '+time_obs,/yoh) ;======================================================================= ;FIRST, find closest MDI files, B_los preferable. ;======================================================================= ; Set 'name' for gif/html logging purposes name=oblk.filename if (keyword_set(suffix)) then begin name = name + suffix endif else begin name=name+'_qnd' endelse dateStr=strmid(oblk.filename,5,13) year=strmid(dateStr,0,4) htmlFile=logFileDir+year+'_pnting_log.html' flagFile=logFileDir+year+'_pnting_FLAGS.html' ;Use date to figure out corresponding MDI orbit number;;;;;;;; mdiclose=mdi_time2file(oblk.date_obs) longbe=strlen(mdiclose) wherebe=strpos(mdiclose,'01d.',/reverse_search) mdiday=strmid(mdiclose,wherebe+4,4) ; wherebe + 4 because starts w/ "0" of 01d, 4 chars long IF (vb ge 2) THEN print, 'MDI DAY corresponding to '+dateStr + ': '+ mdiday IF mdiday eq 0 THEN BEGIN err='date could not be translated to an MDI orbit number' return ENDIF mdipath='/data/SOLDAT/MDI/'+strtrim(year,2)+'/' mdiIdir=mdipath+'fd_Ic_6h_01d.'+mdiday+'/' mdiBdir=mdipath+'fd_M_96m_01d.00'+mdiday+'/' wherebe=strpos(mdiclose,'/',/reverse_search) bestBfile=strmid(mdiclose,wherebe+1,longbe-wherebe) ; get just the filename bestfile=mdiBdir+bestBfile ; full path +filename check=file_search(bestfile) ; make sure that gram actually exists IF (strlen(check) lt 2) THEN BEGIN ; if file doesn't exist, use closest 'gram on same day. anygram=file_search(mdiBdir+'/fd_M_96m_01d*.0*.fits') IF strlen(anygram(0)) gt 0 THEN BEGIN whichwas = fix(strmid(bestbfile,18,4)) ;whichanygram=where(anygram eq bestfile) closeorder=abs(indgen(n_elements(anygram))-whichwas) closeorder=sort(closeorder) bestfile=anygram(closeorder(0)) ENDIF ELSE BEGIN ; then try for an Ic file ; logic: use Ic file 0001 if closest B-file was in first half of the day, ; Ic file 0003 if in second half. keep it simple use_blong = 0 whichIc='0000' whichbe=fix(strmid(bestBfile,8,4,/rev)) ; just the 0013, or similar, file number. dayfrac=whichbe/15. ; there are 15 B-grams per MDI day IF (dayfrac ge 0.20) then whichIc='0001' IF (dayfrac ge 0.45) then whichIc='0002' IF (dayfrac ge 0.70) then whichIc='0003' bestIfile='fd_Ic_6h_01d.'+mdiday+'.'+whichIc+'.fits' bestfile=mdiIdir+bestIfile check=file_search(bestfile) IF (strlen(check) lt 2) THEN BEGIN ; if that nice way works, use any continuum file found. anyicont=file_search(mdiIdir+'/fd_Ic*_01d.*.0*.fits') IF strlen(anyicont(0)) gt 0 THEN BEGIN bestfile=anyicont(0) ENDIF ELSE BEGIN print,'WARNING, NO MDI DATA' err = 'WARNING, NO MDI DATA' return ENDELSE ENDIF ENDELSE ENDIF ;======================================================================= ; Get MDI Data, and correct for rotation ;======================================================================= if (vb ge 2) then quiet=0 else quiet=1 mreadfits, bestfile, mindex, mdata, quiet=quiet nan_pos = where( finite( mdata ) eq 0, nnans) badfrac=nnans/float(n_elements(mdata)) IF badfrac gt 0.5 THEN BEGIN anygram=file_search(mdiBdir+'/fd_M_96m_01d*.0*.fits') whichanygram=where(anygram eq bestfile) closeorder=abs(indgen(n_elements(anygram))-whichanygram[0]) closeorder=sort(closeorder) i=1 WHILE ((badfrac gt 0.5) && (i lt n_elements(anygram)-1)) DO BEGIN mreadfits, anygram(closeorder(i)), mindex, mdata, quiet=quiet nan_pos = where( finite( mdata ) eq 0 , nnans) badfrac=nnans/float(n_elements(mdata)) i++ ; ELW20131004 ENDWHILE if (badfrac gt 0.5) then begin ; ELW20131004 print,'WARNING, MDI DATA CONTAINS TOO MANY NaNs' err = 'WARNING, MDI DATA CONTAINS TOO MANY NaNs' return endif ENDIF ;correct for MDI's spacecraft yaw/roll angle nxmdi = n_elements(mdata[*,0]) nymdi = n_elements(mdata[0,*]) mdidate = anytim(mindex.date_obs,/yoh) ddays = int2secarr(mdidate,ivmdate0)/86400. IF (abs(mindex.solar_p) gt 40) then BEGIN mindex.solar_p0=mindex.solar_p0 + 180. mindex.solar_p=mindex.solar_p + 180. mindex.p_angle=mindex.p_angle + 180. im=mdata[*,*] im=rotate(im,2) mdata[*,*]=im ENDIF mdipix = mindex.im_scale mdipix_earth = mdipix*mindex.EARTH_R0/mindex.OBS_R0 mdilatlon=get_latlon(nx=nxmdi,ny=nymdi,x0=mindex.xcen,y0=mindex.ycen,$ ; x0,y0 should ; be close to 0,0 b0=mindex.obs_b0*!dtor,pang=mindex.solar_p*!dtor,pixs=mdipix,rad=mindex.obs_r0) mdilatlon.cmd=mdilatlon.cmd+diff_rot(ddays,mdilatlon.latitude,/allen) ;======================================================================= ; Get IVM Data, and correct as needed ;======================================================================= cmd = oblk.clng lat = oblk.clat radius = oblk.radius b0 = oblk.b0 p = oblk.p pix_size = [ oblk.cdelt1, oblk.cdelt2 ] if ~keyword_set(icont) then begin print, 'ivm_mdi_align: Icont not set. Exiting.' return endif nxivm = n_elements(blong[*,0]) nyivm = n_elements(blong[0,*]) xyivm = lonlat2xy([cmd*!radeg,lat*!radeg], b0=b0,radius=radius) xyivmsub='' IF use_blong then BEGIN ivmdat=blong ENDIF ELSE BEGIN if ~keyword_set(icont) then begin print, 'ivm_mdi_align: icont keyword not set. Returning.' return endif ivmdat=icont ENDELSE ;find limb: ;R = IVM_MK_RADIUS( oblk ) ;rmax = R - oblk.radius ; ELW20100325 - closer to orig ; rmax = R - radius ;mnex=-20. ;limb = WHERE( rmax GT mnex, nlimb ) ivm_skymask=fltarr(nxivm,nyivm) xivm=fltarr(2) yivm=fltarr(2) ;trimmed_ivm=ivmdat[trim:nxivm-trim-1, trim:nyivm-trim-1] trimmed_ivm=ivmdat sstrim=size(trimmed_ivm) xivm[0]=0 xivm[1]=sstrim[1]-1 yivm[0]=0 yivm[1]=sstrim[2]-1 ;ivm_skymask=ivm_skymask[trim:(nxivm-trim-1), trim:(nyivm-trim-1)] ;matches ivm images w/ trimmed edges ;========================================================================= ; LOOK FOR CONTRAST IN MDI IMAGE. THEN CROP PART OF IMAGE TO USE FOR ; ALIGNMENT. (INPUT TO CORREL_OPTIMIZE) ;========================================================================= ;input trimmed ivm image ; spotArray=ivm_lookForSunspots(trimmed_ivm, oblk, ivm_skymask, imageDir, name, blong=use_blong, /no_tvrd) ;IF strtrim(spotArray[0],2) eq 'noContrast' then BEGIN ;no contrast in first image checked ;(typically, blong. Try cont. image next) ; print, "WARNING! No significant contrast for alignment! " ; IF use_blong then BEGIN ; im_t = bytscl(trimmed_ivm, -.10, .10) ; ENDIF ELSE BEGIN ; im_t = bytscl(trimmed_ivm, 5000, max(trimmed_ivm)) ; ENDELSE ;ENDIF ELSE BEGIN ;trim skymask to match spotarray rectangle ;ivm_skymask=ivm_skymask[spotArray[0]+trim:spotArray[1]+trim, spotArray[2]+trim:spotArray[3]+trim] ; ivm_skymask=ivm_skymask[spotArray[0]:spotArray[1], spotArray[2]:spotArray[3]] ; xivm[0]=spotArray[0] ; xivm[1]=spotArray[1] ; yivm[0]=spotArray[2] ; yivm[1]=spotArray[3] ;ENDELSE ; Now get the equivalent region in the MDI data ;retrieve same MDI data (scaled to match pixels of IVM) but add ;pixels around all sides so that the MDI image is larger. This is ;necessary for the cross-correlation function to work properly. edge=150 xbigIvm=fltarr(2) ybigIvm=fltarr(2) xbigIvm[0]=xivm[0]-edge xbigIvm[1]=xivm[1]+edge ybigIvm[0]=yivm[0]-edge ybigIVM[1]=yivm[1]+edge ; account for shift in center when using restricted FOV xpshift = (xbigIVM[0]/2.0 + (xbigIvm[1]-(nxivm-1L))/2.0)*pix_size[0] ypshift = (ybigIVM[0]/2.0 + (ybigIvm[1]-(nyivm-1L))/2.0)*pix_size[1] ; convert the shift from image coords to solar coords xshift = xpshift*cos(-p)-ypshift*sin(-p) yshift = xpshift*sin(-p)+ypshift*cos(-p) ;xcen,ycen are the coordinates of the center of the IVM's FOV (arcsec), from old pointing xcen = xyivm[0] + xshift ycen = xyivm[1] + yshift ; Correct these for solar rotation ll = xy2lonlat([xcen,ycen],b0=b0,radius=radius) ll[0] = ll[0] + diff_rot(ddays,ll[1],/allen) xy = lonlat2xy(ll,mdidate) xcen = xy[0] ycen = xy[1] ;stop ; Now convert to MDI pixel coordinates; ; xvi/vyi are LL,UR coords ;XX KDL change to mdipix_earth??? xmdi = fltarr(2) & ymdi = fltarr(2) ;xmdi[0] = mindex.crpix1+xcen/mdipix - $ ;((xbigIvm[1]-xbigIvm[0])/2)*pix_size[0]/mdipix ;xmdi[1] = mindex.crpix1+xcen/mdipix + $ ;((xbigIvm[1]-xbigIvm[0])/2)*pix_size[0]/mdipix ;ymdi[0] = mindex.crpix2+ycen/mdipix - $ ;((ybigIvm[1]-ybigIvm[0])/2)*pix_size[1]/mdipix ;ymdi[1] = mindex.crpix2+ycen/mdipix + $ ;((ybigIvm[1]-ybigIvm[0])/2)*pix_size[1]/mdipix xmdi[0] = mindex.crpix1+xcen/mdipix_earth - $ ((xbigIvm[1]-xbigIvm[0])/2)*pix_size[0]/mdipix_earth xmdi[1] = mindex.crpix1+xcen/mdipix_earth + $ ((xbigIvm[1]-xbigIvm[0])/2)*pix_size[0]/mdipix_earth ymdi[0] = mindex.crpix2+ycen/mdipix_earth - $ ((ybigIvm[1]-ybigIvm[0])/2)*pix_size[1]/mdipix_earth ymdi[1] = mindex.crpix2+ycen/mdipix_earth + $ ((ybigIvm[1]-ybigIvm[0])/2)*pix_size[1]/mdipix_earth ; Shift by the fractional part of xvi,yvi xshift = xmdi[0]-long(xmdi[0]) & yshift = ymdi[0]-long(ymdi[0]) xmdi = long(xmdi) ymdi = long(ymdi) ; Make sure the boundaries are not out of range. IF xmdi[0] LT 0 then BEGIN xshift = xshift - xmdi[0] xmdi[0] = 0 ENDIF IF ymdi[0] LT 0 then BEGIN yshift = yshift - ymdi[0] ymdi[0] = 0 ENDIF IF xmdi[1] GT nxmdi-1 then BEGIN xmdi[1] = nxmdi-1 ENDIF IF ymdi[1] GT nymdi-1 then BEGIN ymdi[1] = nymdi-1 ENDIF ; Get the transformation putting mdi onto ivm. tmdi2ivmguess = rss2pq(xshift=xshift, $ yshift=yshift, $ rot12=p*!radeg-mindex.solar_p, $ ; xscale=mdipix/pix_size[0], $ ; yscale=mdipix/pix_size[1], $ xscale=mdipix_earth/pix_size[0], $ yscale=mdipix_earth/pix_size[1], $ xmdi[1]-xmdi[0]+1,ymdi[1]-ymdi[0]+1,/center) ivmdata=trimmed_ivm[xivm[0]:xivm[1],yivm[0]:yivm[1]] mdidata=poly_2d(mdata[xmdi[0]:xmdi[1],ymdi[0]:ymdi[1]], $ tmdi2ivmguess[*,0],tmdi2ivmguess[*,1],1, $ (xbigIvm[1]-xbigIvm[0]+1),(ybigIvm[1]-ybigIvm[0]+1)) mdilat=poly_2d(mdilatlon.latitude[xmdi[0]:xmdi[1],ymdi[0]:ymdi[1]], $ tmdi2ivmguess[*,0],tmdi2ivmguess[*,1],1, $ (xbigIvm[1]-xbigIvm[0]+1),(ybigIvm[1]-ybigIvm[0]+1)) mdicmd=poly_2d(mdilatlon.cmd[xmdi[0]:xmdi[1],ymdi[0]:ymdi[1]], $ tmdi2ivmguess[*,0],tmdi2ivmguess[*,1],1, $ (xbigIvm[1]-xbigIvm[0]+1),(ybigIvm[1]-ybigIvm[0]+1)) s_ivm=size(ivmdata) s_mdi=size(mdidata) ;define offset of image centers (now that center pix should be same ;lat/lon for each): ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;; ELW-20100422 *************** ;; Should this not factor in the difference between the cropped image and original image center? ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ivm_center=[s_ivm[1]/2., s_ivm[2]/2.] mdi_center=[s_mdi[1]/2., s_mdi[2]/2.] ; ELW-20100428 - Factor in post-crop offset from original reference pixel if keyword_set(ref_offset) then begin x_cr = ref_offset[0] & y_cr = ref_offset[1] endif else begin x_cr = 0 & y_cr = 0 endelse xoffInit=mdi_center[0]-ivm_center[0]+x_cr ; ELW-20100423 Was mdi_center[1] yoffInit=mdi_center[1]-ivm_center[1]+y_cr ;stop ;define mdi skymask ; WARNING - assumes mdi limb data set to NaN, as in for lev1.8 mdi_skymask=fltarr(s_mdi[1], s_mdi[2]) nan_pos = where( finite( mdidata ) eq 0 ) IF nan_pos[0] ne -1 then BEGIN mdi_skymask[nan_pos]=1 mdidata[where(mdi_skymask eq 1)]=0 ;set limb equal to 0 (ivmlimb should already be 0) ENDIF ;set up for cross-correlation function: im1=ivmdata im2=mdidata minsize=min([s_ivm[1], s_ivm[2]]) IF use_blong then BEGIN ; BLONG IMAGES: ; Scale IVM two different ways [using KD's idea] - ELW-20100623 good=where(mask eq 1, count) ; Compute standard deviation im1g=im1[good] ; only on good pixels if (count le 3) then im1g=im1 m1=moment(abs(im1g), sdev=sd1) im1b=bytscl(im1, -10.*sd1, 10.*sd1) im3b=bytscl(im1, -2.*sd1, 2.*sd1) ; im1b=bytscl(im1, -.02, .02) ; ELW20100513.1733 ; im3b=bytscl(im1, -.10, .10) ; ELW20100513.1733 ; Scale MDI two different ways im2b=bytscl(im2, -100, 100) im4b=bytscl(im2, -500, 500) ENDIF ELSE BEGIN ;if comparing CONTINUUM IMAGES: im2=bytscl(im2, 5000, max(im2)) im1=bytscl(im1, 5000, max(im1)) im2b=im2-smoothe(im2,50) IF minsize lt 50 then BEGIN ;sharpen the image ; Made this smoothe smaller b/c images may sometimes be a box with minsize ~50x50 pix. im1b=im1-smoothe(im1, minsize) ENDIF ELSE BEGIN im1b=im1-smoothe(im1, 50) ENDELSE IF n_elements(where(mdi_skymask eq 1)) ne 1 then BEGIN im2b[where(mdi_skymask eq 1)]=0 ;set sky equal to 0 ENDIF IF n_elements(where(ivm_skymask eq 1)) ne 1 then BEGIN im1b[where(ivm_skymask eq 1)]=0 ;set sky equal to 0 ENDIF ENDELSE buffer=10 ;correl_optimize, im1b,im2b, x, y, numpix=sqrt(sqrt(s_mdi[0]*s_mdi[1])) correl_optimize, im1b,im2b, x, y, /numpix ; ELW-20100415 ;compute shift in images (cross-correl function) IF (use_blong) then BEGIN ; correl_optimize, im3b,im4b, xcheck, ycheck, numpix=sqrt(sqrt(s_mdi[0]*s_mdi[1])) correl_optimize, im3b,im4b, xcheck, ycheck, /numpix ; ELW-20100415 if (vb ge 2) then print, '**** x=',x, ' xcheck=',xcheck, ' y=',y, ' ycheck=',ycheck IF xcheck gt (x+buffer) || xcheck lt (x-buffer) || ycheck gt (y+buffer) || ycheck lt (y-buffer) then BEGIN print, 'WARNING! Shifts different if bytscled differently. May be failure.' IF x gt 0 || y gt 0 && xcheck lt 0 && ycheck lt 0 then BEGIN x=xcheck ;x & y should always both be negative at ;this point since im2 is bigger than im1. y=ycheck ENDIF ;ELW-20100427 Changed this from a print warning to an error err = err+'Probable failure. Shifts different if bytscled differently.' pntng_chng=0 if keyword_set(noflag) then return ; ELW-20100427 ENDIF ENDIF IF x gt 0 or y gt 0 then BEGIN err=err+"Failure! Shift detected in wrong direction." pntng_chng=0 if keyword_set(noflag) then return ; ELW-20100427 ENDIF ; NOTE: x and y are the shifts to apply to mdidata to line it up with ivmdata, ; where ivmdata is a cropped/smaller array as mdidata. We want to know, how ; big would this shift be if the images had been identical in size to BEGIN ; with. (Eliminate part of shift due to image size difference). ;if images were the same size, real offset would be: xTrue=xoffInit+x ;xoffinit & yoffInit are the pixel offsets of the 2 image centers where centers are same lat/lon. yTrue=yoffInit+y ;IF abs(xTrue) gt 200 OR abs(yTrue) gt 300 then BEGIN IF abs(xTrue) ge 250 OR abs(yTrue) gt 300 then BEGIN ; ELW-20100426 err=err+"Shift unusually large. Check! Pointing NOT changed in file." pntng_chng=0 if keyword_set(noflag) then return ; ELW-20100427 ENDIF ;ivmF, mdiF are images that are the same size that have the same center in terms of lat/lon: ivmF=ivmdata shifted_mdiF=fshift(mdidata, x, y) shifted_mdiF=shifted_mdiF[0:(s_mdi[1]-1)-xoffInit, 0:(s_mdi[2]-1)-yoffInit] shifted_mdilat=fshift(mdilat, x, y) shifted_mdilat=shifted_mdilat[0:(s_mdi[1]-1)-xoffInit, 0:(s_mdi[2]-1)-yoffInit] shifted_mdicmd=fshift(mdicmd, x, y) shifted_mdicmd=shifted_mdicmd[0:(s_mdi[1]-1)-xoffInit, 0:(s_mdi[2]-1)-yoffInit] ;by definition, xtrue, ytrue are the shifts performed on mdiF to match ;ivmF. So, shift on ivmF to match mdiF is then -xtrue, -ytrue ;CALCULATE NEW POINTING: oldlat=lat oldcmd=cmd oldxy=hel2arcmin(lat/!dtor,cmd/!dtor,date=oblk.date_obs) * 60. IF pntng_chng then BEGIN newlat=shifted_mdilat(ivm_center[0],ivm_center[1])*!dtor newcmd=shifted_mdicmd(ivm_center[0],ivm_center[1])*!dtor ; newxy=hel2arcmin(newlat/!dtor,newcmd/!dtor,date=oblk.date_obs,p=0.,b0=0.)*60. txx=cos(newlat)*sin(newcmd)*oblk.radius tyy=(sin(newlat)*cos(oblk.b0)-cos(newlat)*cos(newcmd)*sin(oblk.b0))*oblk.radius txcen=(txx*cos(oblk.p)-tyy*sin(oblk.p)) tycen=(txx*sin(oblk.p)+tyy*cos(oblk.p)) newxy=[txcen,tycen] ; doing p=0, b0=0 seems screwy but it's the magic of *really* going back to IVM frame of ref ; from MDI. ENDIF ELSE BEGIN newlat=oldlat newcmd=oldcmd newxy=oldxy ENDELSE IF (vb ge 3) && pntng_chng THEN BEGIN window,3,xsize=512,ysize=512 loadct,0 tv,rebin(bytscl(rot(mdata,-1*oblk.p/!dtor),-500,500),512,512) loadct,5 draw_grid,[0,0],[2.1,2.1],latlon=[oldlat/!dtor,oldcmd/!dtor],$ tilt=oblk.b0/!dtor,roll=oblk.p/!dtor,gridsep=10.,col=100,/noer draw_grid,[0,0],[2.1,2.1],latlon=[newlat/!dtor,newcmd/!dtor],$ tilt=oblk.b0/!dtor,roll=oblk.p/!dtor,gridsep=10.,/noer,col=200 loadct,0 ; tv,bytscl(ivmF,-0.05,0.05),3 window,4,xsize=(size(ivmF))[1],ysize=(size(ivmF))[2] tv,bytscl(ivmF,-0.15,0.15) wait, 2 ;stop wdelete,3 wdelete,4 ENDIF ; window,4,xsize=850,ysize=512 ; ELW-20100312 if ((vb ge 3) || keyword_set(log)) then begin shifts=[xTrue, yTrue] files=[htmlFile, flagFile] IF use_blong then BEGIN if (vb ge 1) then print, 'Calling ivm_savegifs (use_blong)...' saved=ivm_savegifs(imageDir, files, name, ivmF, shifted_mdiF, shifts, $ oldlat, oldcmd, [newlat, newcmd], warning=err, /blong, /qnd, /eps) ENDIF ELSE BEGIN if (vb ge 1) then print, 'Calling ivm_savegifs (default)...' saved=ivm_savegifs(imageDir, files, name, ivmF, shifted_mdiF, shifts, $ oldlat, oldcmd, [newlat, newcmd], warning=err) ENDELSE endif if (vb ge 2) then begin IF (oldlat eq newlat) && (newcmd eq oldcmd) then BEGIN print, "WARNING: No change made to pointing!" ENDIF ELSE BEGIN if (vb ge 3) then begin s=size(shifted_mdiF) erase IF use_blong then BEGIN tvscl, bytscl(shifted_mdiF, -500, 500), 0, 0 tvscl, bytscl(ivmF,-500,500), 0,s[2] ENDIF ELSE BEGIN tvscl, shifted_mdiF, 0, 0 tvscl, ivmF, 0,s[2] ENDELSE endif ;ELW20110218 XXXX Check ddays, output newlat, newcmd too to see if something weird is going on print, date_obs,time_obs,' ',strtrim(ddays,2),' ',bestBfile, newlat,newcmd print, 'oldlat, oldcmd,newlat,newcmd, newxy = ', oldlat,oldcmd,newlat, newcmd, newxy print, 'Old center: N/S '+strtrim(oldlat/!dtor,2)+' E/W '+strtrim(oldcmd/!dtor,2) print, 'New center: N/S '+strtrim(newlat/!dtor,2)+' E/W '+strtrim(newcmd/!dtor,2) print, 'shift: '+strtrim(sqrt((newxy[0]-oldxy[0])^2+(newxy[1]-oldxy[1])^2),2) +' arcsec' ENDELSE endif END