; docformat = 'rst' ; ; NAME: ; showAoImages ; PURPOSE: ; Create a figure of multiple AO images ;+ ; :Description: ; Create a figure of multiple AO images ; ; :Categories: ; Figures ; ; :Params: ; file: in, required, type=string ; Path to the file listing plot parameters ; ; :Keywords: ; ; :Uses: ; saiLoadConfig ; ; :Author: ; B.Carry (OCA) ; ; :History: ; Change History:: ; Original Version written in 2016, B. Carry (OCA):: ;- function showAoImages, file COMPILE_OPT hidden, idl2 ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- I -- Initialization And Input Verification -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--I.1-- Test Existance of Input File -------------------------------------------------------- if not keyword_set(file) then $ if file_test('plot.ini',/Read) then file='plot.ini' else begin message, /ioError, 'Cannot find the input file: '+strtrim(file,2) return, -1 ;--Cannot read input file endelse ;--I.2-- Read Configuration File ------------------------------------------------------------- conf = saiLoadConfig( file ) sel = where( conf.image.show eq 1, nbSel ) ;--I.X-- TBD --------------------------------------------------------------------------------- facZ = 8 facPS= 500 key={xLeft:0.05, xRight:0.95, $ yTop:0.95, yShi:0.075, $ comp:{x:0.85, y:0.10}, $ bar:{x:0.25, y:0.10}} theta = 2*!PI*fIndGen(361)/360 charSize = 4 charThick= 7 charCol ='White' ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- II -- Definition of Graphics Output Parameters -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--II.1-- Save Current Device Status --------------------------------------------------------- SetDecomposedState, 1, CurrentState=theState thisDevice = !D.Name ;--II.2-- Layout Parameters for Figures ------------------------------------------------------ ;--II.2.1-- Plot Positions xBot = 50. ;-- Left Margin yBot = 50. ;-- Bottom Margin xLength = 8000. ;-- Plot X Size yLength = 8000. ;-- Plot Y Size xTop = 50. ;-- Right Margin yTop = 50. ;-- Top Margin xSep = 50 ;-- Plot-to-plot X Separation ySep = 50 ;-- Plot-to-plot Y Separation ;--II.2.2-- Canvas Plot Size winZ = [ xBot + xTop + (conf.graph.N.col-1)*(xSep+xLength) + xLength, $ yBot + yTop + (conf.graph.N.row-1)*(ySep+yLength) + yLength ] /facZ winPS = winZ *facZ / facPS ;--II.3-- Loop on Selected Images: Open/Close Figures --------------------------------------- winInd = indgen( ceil(nbSel / float(conf.graph.N.row*conf.graph.N.col)) )+1 nbPlotted = 0 for kSel=0, nbSel-1 do begin print, 'Image '+strTrim(string(kSel+1),2)+'/'+strTrim(string(nbSel),2) ;--II.3.1-- A new File must be Opened if (float(kSel-nbPlotted) mod (conf.graph.N.row*conf.graph.N.col)) eq 0 or $ kSel eq 0 then begin ;--II.3.2-- Close Former Plot if kSel gt 0 then begin snapshot = cgSnapshot() cgPS_open, fileName=conf.graph.name+'-'+string(winInd[kWin],format='(I03)')+'.eps', $ /Metric, /Decomposed, /Portrait, $ xSize=winPS[0], ySize=winPS[1], language_level=2, /quiet cgImage, snapshot cgPS_close, /png, Delete_PS=0 endif ;--II.3.3-- Update the Number of Columns and Rows nbPlotted = kSel kWin = kSel/(conf.graph.N.row*conf.graph.N.col) trialRow = ceil( (nbSel-kSel) / conf.graph.N.col)+1 trialCol = ceil( (nbSel-kSel) mod conf.graph.N.col) if trialCol eq 0 then trialRow-- if trialRow lt conf.graph.N.row then conf.graph.N.row=trialRow if conf.graph.N.row eq 1 then begin if trialCol lt conf.graph.N.col and trialCol ne 0 then conf.graph.N.col=trialCol endif ;--II.3.4-- Update the Figure Size winZ = [ xBot + xTop + (conf.graph.N.col-1)*(xSep+xLength) + xLength, $ yBot + yTop + (conf.graph.N.row-1)*(ySep+yLength) + yLength ]/facZ winPS = winZ *facZ / facPS ;--II.3.5-- Open a New Z-buffer for current views Set_Plot, 'Z' device, Set_Resolution=winZ, Decomposed=1, Set_Pixel_Depth=24 cgErase endif ;--II.4-- Index of Current Image in Grid --------------------------------------------------- kCol = ((kSel-nbPlotted) mod (conf.graph.N.row*conf.graph.N.col)) mod conf.graph.N.col kRow = ((kSel-nbPlotted) mod (conf.graph.N.row*conf.graph.N.col)) / conf.graph.N.col ;--II.5-- Select Current Image ------------------------------------------------------------- loc = conf.image[sel[kSel]] ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- III -- Read Image and Extract Parameters -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--III.1-- Read the Image ------------------------------------------------------------------ image = readFits( loc.name, header, /Silent ) dimI = size(image) ;--III.2-- Get Observatory, Camera, etc ---------------------------------------------------- tel = getTelescopeId(header) cam = getCameraId(header) filter = getCameraFilter(header, tel, cam ) pixel = getCameraPixel( header, tel, cam ) pixel=pixel.resol ;--III.3-- Get Observing Conditions: Epoch, Orientation, etc ------------------------------- dit = getCameraDIT( header, tel, cam ) pa = getCameraPA( header, tel, cam ) ;--III.4-- Define Window Sizes ------------------------------------------------------------- ;--III.4.1-- Main Window case loc.main.unit of 'pix': mainSize = ceil(loc.main.size) 'sec': mainSize = ceil(loc.main.size/pixel) 'min': mainSize = ceil(loc.main.size/pixel) /60. 'deg': mainSize = ceil(loc.main.size/pixel) /3600. ; 'km' : stop ;-TBD endcase ;--III.4.2-- Inset Window case loc.inset.unit of 'pix': insetSize = ceil(loc.inset.size) 'sec': insetSize = ceil(loc.inset.size/pixel) 'min': insetSize = ceil(loc.inset.size/pixel) /60. 'deg': insetSize = ceil(loc.inset.size/pixel) /3600. endcase ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- IV -- Display the Main Image -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--IV.1-- Valid Layer ---------------------------------------------------------------------- if (loc.main.layer ne 0 and dimI[0] ne 3) or $ (loc.main.layer gt dimI[dimI[0]]-1) then begin message, /Informational, 'Invalid layer (#'+strTrim(string(loc.main.layer),2)+') for image '+loc.name+' ('+loc.nick+')' set_plot, thisDevice SetDecomposedState, theState return, -11 ;--Invalid layer for main image endif ;--IV.2-- Extract and Rotate the View ------------------------------------------------------ if strCmp(loc.main.track,'det') then begin main=cropFrame( image[*,*,loc.main.layer], [dimI[1],dimI[2]]/2., mainSize ) print, 'main no rotation!' endif else begin main=cropFrame( image[*,*,loc.main.layer], [dimI[1],dimI[2]]/2., mainSize, angle=pa, header=header) endelse ;--IV.3-- Location of Current View --------------------------------------------------------- posView = [xBot + kCol*(xSep+xLength), $ yBot + (conf.graph.N.row-1-kRow)*(ySep+yLength), $ xBot + kCol*(xSep+xLength) + xLength, $ yBot + (conf.graph.N.row-1-kRow)*(ySep+yLength) + yLength ] posView[[0,2]] /= (winZ[0]*facZ) posView[[1,3]] /= (winZ[1]*facZ) ;--IV.4-- Display Current View ----------------------------------------------------------- cgImage, main, /NoErase, position=posView, $ minValue=loc.main.min, maxValue=loc.main.max, CTIndex=loc.main.CTIndex ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- V -- Display an Optional Inset -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; if loc.inset.layer ne -1 then begin ;--V.1-- Valid Layer --------------------------------------------------------------------- if (loc.inset.layer ne 0 and dimI[0] ne 3) or $ (loc.inset.layer gt dimI[dimI[0]]-1) then begin message, /Informational, 'Invalid layer (#'+strTrim(string(loc.inset.layer),2)+') for image '+$ loc.name+' ('+loc.nick+')' set_plot, thisDevice SetDecomposedState, theState return, -12 ;--Invalid layer for inset image endif ;--V.2-- Extract and Rotate the View ----------------------------------------------------- if strCmp(loc.main.track,'det') then begin inset = cropFrame( image[*,*,loc.inset.layer], [dimI[1],dimI[2]]/2., mainSize ) print, 'inner disk no rotation!' endif else begin inset = cropFrame( image[*,*,loc.inset.layer], [dimI[1],dimI[2]]/2., mainSize, angle=pa ) endelse ;--V.3-- Define the Inset Region of Interest --------------------------------------------- distArr=shift( dist( mainSize, mainSize ), mainSize/2., mainSize/2. ) roi = where( distArr le insetSize, comp=mask ) ;--V.4-- Create a Transparent Image ------------------------------------------------------ miss = 10101 inset[mask] = miss ;--V.5-- Display Current View ------------------------------------------------------------ cgImage, inset, /NoErase, alphaFGPosition=posView, $ CTIndex=loc.inset.CTIndex, Transparent=0, Missing_Value=miss endif ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- VI -- Display of Satellites in Multiple Systems -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; if loc.label.sat eq 1 then begin ;--VI.1-- Get Satellite Ephemerides ------------------------------------------------------; tNow=string(systime(/UTC,/JULIAN),format='(D16.8)') tmp1='/tmp/'+tNow+'-1.dat' tmp2='/tmp/'+tNow+'-2.dat' ;--VI.1.1-- First Satellite print, loc.satellite.S1 if not strCmp(loc.satellite.S1,'') then $ eSat1 = voMiriade_callEphemCC( loc.satellite.S1, dit.start, /Aster, /Web, /ISO, dump=tmp1 ) ;--VI.1.2-- Second Satellite if not strCmp(loc.satellite.S2,'') then $ eSat2 = voMiriade_callEphemCC( loc.satellite.S2, dit.start, /Aster, /Web, /ISO, dump=tmp2 ) print, dit.start if strcmp(loc.satellite.s1,'379/1') then begin print, 'ATTENTION!!!! Teak momentatee pour HUENNA. a virer!!!!' case dit.start of '2004-12-10T06:51:28.80': begin print, 'VLT2004' print, eSat1[0].dx, eSat1[0].dy eSat1[0].dx = 1.650 eSat1[0].dy = 0.160 end '2005-01-26T02:47:39.26': begin print, 'VLT2005 - Jan' print, eSat1[0].dx, eSat1[0].dy eSat1[0].dx = -1.888 eSat1[0].dy = -0.09 end '2005-02-02T03:09:16.21': begin print, 'VLT2005 - Feb' print, eSat1[0].dx, eSat1[0].dy end '2005-02-16T01:21:08.59': begin print, 'VLT2005 - Feb16' print, eSat1[0].dx, eSat1[0].dy eSat1[0].dx = 1.266 eSat1[0].dy = 0.00 end '2009-11-28T12:13:35.10': begin print, 'keck 2009 11' print, eSat1[0].dx, eSat1[0].dy end '2012-02-03T13:18:43.99': begin print, 'keck 2012 feb' print, eSat1[0].dx, eSat1[0].dy end '2012-03-04T10:33:13.88': begin print, 'keck 2012 mar' print, eSat1[0].dx, eSat1[0].dy eSat1[0].dx = -1.200 eSat1[0].dy = 0.894 print, eSat1[0].dx, eSat1[0].dy end else: endcase endif ;--VI.2-- Set Up Plane of the Sky Projection ---------------------------------------------; case loc.main.unit of 'pix': fovSize = 0.5* loc.main.size*pixel 'sec': fovSize = 0.5* loc.main.size 'min': fovSize = 0.5* loc.main.size*pixel *60. 'deg': fovSize = 0.5* loc.main.size*pixel *3600. endcase ;--VI.3-- Plot the Satellite Positions ---------------------------------------------------; ;--VI.1.1-- First Satellite if not strCmp(loc.satellite.S1,'') then begin cgPlot, eSat1[0].dx, eSat1[0].dy, /NoErase, position=posView, $ xStyle=1+4, xRange=fovSize*[ 1,-1], $ yStyle=1+4, yRange=fovSize*[-1,1 ], $ psym='Open Circle', symSize=10, thick=6 endif ;--VI.1.2-- Second Satellite if not strCmp(loc.satellite.S2,'') then begin cgPlot, eSat2[0].dx, eSat2[0].dy, /NoErase, position=posView, $ xStyle=1+4, xRange=fovSize*[ 1,-1], $ yStyle=1+4, yRange=fovSize*[-1,1 ], $ psym='Open Circle', symSize=10, thick=3 endif endif ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- VII -- Display Labels, Compass, Scale bar -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--VI.1-- Set the Graphic Area ------------------------------------------------------------- cgPlot, 0,0, /NoErase, /NoData, position=posView, $ xStyle=1+4, xRange=[0,1], $ yStyle=1+4, yRange=[0,1] ;--VII2-- Mark Inset Limit if Needed ------------------------------------------------------- if loc.inset.layer ne -1 then begin x = (1.*insetSize/mainSize) * cos(theta) + 0.5 y = (1.*insetSize/mainSize) * sin(theta) + 0.5 cgPlot, /OverPlot, x, y, thick=5;, color='black' endif ;--VII.3-- Facility Labels ------------------------------------------------------------------ ;--VII.3.1-- Telescope if loc.label.tel eq 1 then begin cgText, key.xLeft, key.yTop-0*key.yShi, align=0, tel, charSize=charSize, charThick=charThick, color=charCol endif ;--End of VI.3.1 ;--VII.3.2-- Instrument if loc.label.ins eq 1 then begin cgText, key.xLeft, key.yTop-1*key.yShi, align=0, cam, charSize=charSize, charThick=charThick, color=charCol endif ;--End of VI.3.2 ;--VII.3.3-- Filter if loc.label.ins eq 1 then begin cgText, key.xLeft, key.yTop-2*key.yShi, align=0, filter.id, charSize=charSize, charThick=charThick, color=charCol endif ;--End of VI.3.2 ;--VII.4-- Information on the Observations -------------------------------------------------- ;--VII.4.1-- Date of observation if loc.label.date eq 1 then begin date = strMid( dit.start,0,10) cgText, key.xRight, key.yTop-0*key.yShi, align=1, date, charSize=charSize, charThick=charThick, color=charCol endif ;--End of VII.4.1 ;--VII.4.2-- UT of observation if loc.label.UTC eq 1 then begin ut = strMid( dit.start,11,5) cgText, key.xRight, key.yTop-1*key.yShi, align=1, ut, charSize=charSize, charThick=charThick, color=charCol endif ;--End of VII.4.1 ;--VII.4.3-- Scale Bar if not strCmp(loc.label.bar.unit,'none') then begin case loc.main.unit of 'pix': fovSize = 1.*ceil(loc.main.size)*pixel 'sec': fovSize = 1.*ceil(loc.main.size) endcase case loc.label.bar.unit of 'pix': begin barLen = loc.label.bar.val / fovSize barStr = string(loc.label.bar.val,format='(F3.1)')+' pix' end 'sec': begin barLen = loc.label.bar.val / fovSize barStr = string(loc.label.bar.val,format='(F3.1)')+'"' end endcase cgErrPlot, /Horizontal, key.bar.y, key.bar.x - barLen/2., key.bar.x + barLen/2., $ thick=5, width=0.03, color=charCol cgText, key.bar.x, key.bar.y+0.01, barStr, align=0.5, charSize=charSize, charThick=charThick, color=charCol endif ;--End of VII.4.3 ;--VII.4.4-- Field Orientation if loc.label.sky eq 1 then begin arrows, header, key.comp.x, key.comp.y, /Data, thick=charThick, charSize=charSize, arrowlen=5, color=charCol endif ;--End of VII.4.4 ;--VII.5-- Information on the Target -------------------------------------------------------- if loc.label.name eq 1 then begin cgText, key.xRight, key.yTop-2*key.yShi, align=1, loc.target, charSize=charSize, charThick=charThick, color=charCol endif ;--VII.6-- Panel Id ------------------------------------------------------------------------- if loc.label.panel eq 1 then begin cgText, 0.5, key.yTop-key.yShi/2., strTrim(string(kSel+1),2), align=0.5, $ charSize=charSize, charThick=charThick cgPlot, 0.5, key.yTop-key.yShi/3., /NoErase, position=posView, $ xStyle=1+4, xRange=[0,1], $ yStyle=1+4, yRange=[0,1], $ psym='Open Circle', symSize=10, thick=3 endif endfor ;--II-3--End of Loop on Selected Images ;--II.4-- Close Last Plot -------------------------------------------------------------------- snapshot = cgSnapshot() cgPS_open, fileName=conf.graph.name+'-'+string(winInd[kWin],format='(I03)')+'.eps', $ /Metric, /Decomposed, /Portrait, $ xSize=winPS[0], ySize=winPS[1], language_level=2, /quiet cgImage, snapshot cgPS_close, /png, Delete_PS=0 set_plot, thisDevice SetDecomposedState, theState return, 1 end