; docformat = 'rst' ; ; NAME: ; forcAstronomy ; PURPOSE: ; Long-term planning of astronomical observations ;+ ; :Description: ; Long-term planning of astronomical observations ; ; :Categories: ; Ephemeris, VO, forcAstronomy ; :Params: ; path: in, optional, type=string ; The directory containing forcAstronomy ephemeris file ; ; :Returns: A flag indicative of the completions of tasks:: ; ; :Uses: ; ; :Author: ; B.Carry (OCA) ; ; :History: ; Change History:: ; Original Version written in November 2016, B. Carry (ESA) ;- function forecAstronomy, path, mime=mime, cuts=cuts, conf=conf, debug=debug ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- I -- Initialization And Input Verification -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; COMPILE_OPT hidden, idl2 !Quiet=1 t0 = systime(/JULIAN, /UTC) ;--I.1-- Input Directory --------------------------------------------------------------------- if not keyword_set(path) then path ='./' else path=strTrim(string(path,format='(A)'),2)+'/' if not file_test(path) then begin message, 'Target directory does not exist: '+strtrim(path,2) return, 2 endif ;--I.2-- Input Ephemeris --------------------------------------------------------------------- ;-tbd ephemFile = path+'epheforcast.dat' ;-tbd if not file_test(ephemFile,/read) then begin ;-tbd message, 'Ephemeris file not found: '+strtrim(ephemFile,2) ;-tbd return, 1 ;-tbd endif ;--I.3-- Default Configuration --------------------------------------------------------------- ;--I.3.1-- Forecast Sources Directory forecastPATH = routine_filepath( 'forecastronomy', /Is_Function ) forecastPATH = strMid( forecastPATH, 0, strlen(forecastPATH)-18) ;--I.3.2-- File Extensions sufEPS = '.eps' sufPNG = '.png' ;--I.3.3-- Log location if not keyword_set(debug) then debug='/dev/null' ;--I.4-- Read Configuration File ------------------------------------------------------------- confDefault = forecastLoadConfig(forecastPATH+'config.ini') ;--I.4.1-- Use Default Configuration if not keyword_set(conf) then conf=confDefault else begin typeC=size(conf,/type) ;--I.4.2-- Use Provided Configuration if typeC eq 8 then begin conf = updateStructure( confDefault, conf ) ;--I.4.3-- Read & Use Provided Configuration endif else begin conf = forecastLoadConfig( conf ) conf = updateStructure( confDefault, conf ) endelse endelse ;--I.5-- Selection Criteria ------------------------------------------------------------------ ;--I.5.1-- Complete Default Criteria if not keyword_set(cuts) then cuts = conf.select $ else cuts = updateStructure( conf.select, cuts ) ;--I.6-- Parameter, Display, and Selection Indices ------------------------------------------- ;--I.6.1-- Define Indices for each Parameter ind={vMag:0, el:1, phase:2, elong:{sun:3, moon:4}, dObs:5, dSun:6, $ cEQ:{ra:7,dec:8}, cEc:{lon:9,lat:10}, cGal:{lon:11,lat:12}, tot:13} ;--I.6.2-- Retrieve Set of Parameters to Display dispId=-1 if conf.disp.el eq 1 then dispId=[dispId,ind.el] if conf.disp.vmag eq 1 then dispId=[dispId,ind.vmag] if conf.disp.cEq.dec eq 1 then dispId=[dispId,ind.cEq.dec] if conf.disp.dSun eq 1 then dispId=[dispId,ind.dSun] if conf.disp.dObs eq 1 then dispId=[dispId,ind.dObs] if conf.disp.phase eq 1 then dispId=[dispId,ind.phase] if conf.disp.elong.sun eq 1 then dispId=[dispId,ind.elong.sun] if conf.disp.elong.moon eq 1 then dispId=[dispId,ind.elong.moon] dispId=dispId[1:n_elements(dispId)-1] nbDisp = n_elements(dispId) ;--I.6.x-- Set Selection Parameters selectId=[ind.vMag, ind.dSun, ind.dObs, ind.phase] tI = systime(/JULIAN, /UTC) ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- II -- Read Miriade/ForecAstronomy Ephemeris -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;-tbd inf = 'target.list' readcol, path+'/'+inf, intarget, format='(A)', /Silent nbTarget = n_elements(intarget) nbd= 1825 obs={lon:6., lat:+43.7, name:'Caussols'} eTarg = {type:'', num:0, name:''} target = replicate( eTarg, nbTarget ) for kT=0, nbTarget-1 do begin split=strSplit(intarget[kT],':',/Extract) target[kT].type = split[0] target[kT].name = split[1] endfor eEphem={jd: 0.D, $ ;-Julian date iso: '', $ ;-ISO date vMag: 0., $ ;-Apparent V magnitude diam: 0., $ ;-Apparent diameter cEQ: {ra: 0., $ ;-Right Ascencion dec:0.}, $ ;-Declination cEC: {lon:0., $ ;-Ecliptic longitude lat:0.}, $ ;-Ecliptic latitude cGal: {lon:0., $ ;-Galactic longitude lat:0.}, $ ;-Galactic latitude cHor: {az: 0., $ ;-Local Azimuth alt:0.}, $ ;-Local Altitude motion: 0., $ ;-Apparent motion (rate) dObs: 0., $ ;-Range to observer dSun: 0., $ ;-Heliocentric distance phase: 0., $ ;-Solar Phase Angle elong:{sun: 0., $ ;-Solar Elongation moon:0.} } ;-Lunar Elongation ephem=replicate(eEphem,nbTarget,nbd) for kT=0, nbTarget-1 do begin loc = voMiriade_readEphemcc(path+'/'+intarget[kT]+'.xml', mime='votable', tCoor=6) ephem[kT,*].jd = transpose(loc.jd) ephem[kT,*].iso = transpose(loc.iso) ephem[kT,*].cEQ.ra = transpose(loc.raJ2000.dec) ephem[kT,*].cEQ.dec = transpose(loc.decJ2000.dec) ephem[kT,*].cHor.alt = abs(obs.lat - ephem[kT,*].cEQ.dec) ephem[kT,*].dObs = transpose(loc.dObs) ephem[kT,*].dSun = transpose(loc.dSun) ephem[kT,*].Vmag = transpose(loc.mag) ephem[kT,*].phase = transpose(loc.phase) ephem[kT,*].elong.sun = transpose(loc.elong) ;-conversions for kEp=0, nbd-1 do begin ;--II.7.1-- Equatorial to Ecliptic cEC = frameCoord_Eq2Ec( ephem[kT,kEp].cEQ.ra, ephem[kT,kEp].cEQ.dec ) ephem[kT,kEp].cEC.lon = cEC[0] ;-EC: Longitude ephem[kT,kEp].cEC.lat = cEC[1] ;-EC: Latitude ;--II.7.2-- Equatorial to Galactic cGal = frameCoord_Eq2Gal( ephem[kT,kEp].cEQ.ra, ephem[kT,kEp].cEQ.dec ) ephem[kT,kEp].cGal.lon = cGal[0] ;-Galac: Longitude ephem[kT,kEp].cGal.lat = cGal[1] ;-Galac: Latitude endfor endfor tII = systime(/JULIAN, /UTC) ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- III -- Define Visibility Periods -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--III.1-- Select Time Period based on Obs. Criteria ----------------------------------------- sel = forecastSelect( target, ephem, cuts, selectId ) ;--III.2-- Draw list of Obs. Period ---------------------------------------------------------- run = forecastRun( sel,reform(ephem[0,*].jd) ) ;--III.3-- Output List of Run ---------------------------------------------------------------- for kT=0, nbTarget-1 do begin print, target[kT].type, target[kT].name for kR=0, run[kT].N-1 do begin print, 'Run '+strTrim(string(kR+1),2)+': '+$ (*run[kT].dates)[kR].iso1+' - '+(*run[kT].dates)[kR].iso2 endfor endfor tIII = systime(/JULIAN, /UTC) ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- IV -- Summary Graph Generation -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--IV.1-- Number of Visible Targets ---------------------------------------------------------- nbVis = total(sel[*,*,ind.tot],1) ;--IV.2-- Plot Summary Figure ---------------------------------------------------------------- nameFig = 'forecAstronomy.eps' forecastOpenGraphic, 1, path+nameFig, conf=conf forecastPlotSummary, reform(ephem[0,*].jd), nbVis, conf=conf forecastCloseGraphic tIV = systime(/JULIAN, /UTC) ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- V -- Graphic Generation for each Target -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; for kT=0, nbTarget-1 do begin ;--V.1-- Open Graphic Environment ---------------------------------------------------------- nameFig = target[kt].type+'-'+target[kt].name+'.eps' forecastOpenGraphic, nbDisp, path+nameFig, conf=conf kPlot=0 ;--V.2-- Identify Observability Periods ---------------------------------------------------- obs = where( sel[kT,*,ind.tot] eq 1 ) ;--V.3.1-- Apparent Magnitude -------------------------------------------------------------- if conf.disp.vMag eq 1 then begin loc = where( sel[kT,*,ind.vMag] eq 1 ) forecastPlot, reform(ephem[kT,*].jd), reform(ephem[kT,*].vmag), kPlot, $ selected=loc, observable=obs, conf=conf, $ yLabel='V (mag)' kPlot++ endif ;--V.3.2-- Phase Angle --------------------------------------------------------------------- if conf.disp.phase eq 1 then begin loc = where( sel[kT,*,ind.phase] eq 1 ) forecastPlot, reform(ephem[kT,*].jd), reform(ephem[kT,*].phase), kPlot, $ selected=loc, observable=obs, conf=conf, $ yLabel='Phase angle (!Uo!N)' kPlot++ endif ;--V.3.3-- Peak altitude ------------------------------------------------------------------- if conf.disp.el eq 1 then begin loc = where( sel[kT,*,ind.el] eq 1 ) forecastPlot, reform(ephem[kT,*].jd), reform(ephem[kT,*].cHor.alt), kPlot, $ selected=loc, observable=obs, conf=conf, $ yLabel='Elevation (!Uo!N)' kPlot++ endif ;--V.3.4-- Equatorial Coordinates ---------------------------------------------------------- ;--V.3.4/A-- Declination if conf.disp.cEq.dec eq 1 then begin loc = where( sel[kT,*,ind.cEq.dec] eq 1 ) forecastPlot, reform(ephem[kT,*].jd), reform(ephem[kT,*].cEq.dec), kPlot, $ selected=loc, observable=obs, conf=conf, $ yLabel='Declination (!Uo!N)' kPlot++ endif ;--V.3.5-- Distances ----------------------------------------------------------------------- ;--V.3.5/A-- Range to Observer if conf.disp.dObs eq 1 then begin loc = where( sel[kT,*,ind.dObs] eq 1 ) forecastPlot, reform(ephem[kT,*].jd), reform(ephem[kT,*].dObs), kPlot, $ selected=loc, observable=obs, conf=conf, $ yLabel=cgGreek('Delta')+'!Dobs!N (au)' kPlot++ endif ;--V.3.5/B-- Heliocentric distance if conf.disp.dSun eq 1 then begin loc = where( sel[kT,*,ind.dSun] eq 1 ) forecastPlot, reform(ephem[kT,*].jd), reform(ephem[kT,*].dSun), kPlot, $ selected=loc, observable=obs, conf=conf, $ yLabel=cgGreek('Delta')+'!Dsun!N (au)' kPlot++ endif ;--V.4-- Close Graphic Environment --------------------------------------------------------- forecastCloseGraphic endfor tV = systime(/JULIAN, /UTC) ; label=['Init', 'Option', 'Read', 'Select', 'Summary', 'Graphics'] ; time=[t0,tI,tII,tIII,tIV,tV] ; forprint, label, time, time-t0, (time-t0)*24.*3600, textout=2, format='(A-10,2x,D14.6,2x,F8.6,2x,F7.4)' return, 0 end