; docformat = 'rst' ; ; NAME: ; shapePredictFlux ; PURPOSE: ; Predict the apparent lightcurve of a SSO based on its shape model ;+ ; :Description: ; Predict the apparent lightcurve of a SSO based on its shape model ; ; :Categories: ; Ephemeris, Visualization ; ; :Params: ; ID: in, required, type=integer/string ; Target identifyer in Miriade's format: Name, ; Provisional designation, Number, preceded by "a:" for asteroids, and ; "c:" for comets. ; ep: in, required, type=double/string ; Starting epoch of computation (JD or ISO) or path to a file with multiple epochs. ; shape: in, required, type=string ; Path to a file with the shape model. ; spin: in, required, type=structure ; Path to a file with the spin information. ; ; :Keywords: ; format: in, optional, type=integer, default=1 ; Specify the output format:: ; 1: Simple with magnitude: JD Mag ; 2: Simple with flux: JD Flux ; 3: LCI: JD Flux X_sun Y_sun Z_sun X_obs Y_obs Z_obs ; ; dump: in, optional, type=string ; Path to a file that will contain the lightcurve ; ; nbd: in, optional, type=integer, default=1 ; The requested number of date to compute ; step: in, optional, type=string, default='1h' ; Step of increment followed by one of (d)ays or (h)ours or ; (m)inutes or (s)econds, e.g, '1.5h' ; period: in, optional, type=float ; Target sidereal rotation period (hours). If set, the nbd ; epochs are evenly spread over the period. ; ; observer: in, optional, type=string, default='500' ; The IAU Code or geographic coordinates of the observer's location ; ; :Examples: ; ; :Uses: ; ; ; :Author: ; B.Carry (OCA) ; ; :History: ; Change History:: ; Original Version written in March-May 2016, B. Carry (OCA) ;- function shapePredictFlux, id, ep, shape, spin, format=format, dump=dump, $ nbd=nbd, step=step, period=period, observer=observer ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- I -- Initialization And Input Verification -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; COMPILE_OPT hidden, idl2 ;--I.1-- Input Syntax Verification ----------------------------------------------------------- if N_PARAMS() lt 4 then begin message, /ioError, 'Syntax : LC = shapePredictFlux( id, ep, model, spin [, format=, dump=, nbd=, step=, ...] )' return, -1 endif ;--I.2-- Test if Model and Spin can be Read -------------------------------------------------- if not file_test(shape,/READ) then begin message, /ioError, 'The model cannot be read: '+shape return, -1 endif if not file_test(spin,/READ) then begin message, /ioError, 'The spin cannot be read: '+spin return, -1 endif spinStr = readSpin( spin ) ;--I.3-- Single or Multiple Epochs submitted? ------------------------------------------------ ;if ep = array -> report values for these dates only ;--I.4-- Number of epochs, time step, period ------------------------------------------------- if not keyword_set(nbd) then nbd = 1 if not keyword_set(step) then step = '1h' if keyword_set(period) then begin nbd = 720 step = strtrim(string(spinStr.ph/nbd,format='(F11.6)'),2)+'h' endif ;--I.5-- Output Format ----------------------------------------------------------------------- if not keyword_set(format) then format=1 ;--I.6-- Observing Geometry ------------------------------------------------------------------ if keyword_set(observer) then begin split = strSplit(observer,' ', /EXTRACT, count=N) if N gt 1 then obsName = strTrim(observer,2) $ else obsName = strTrim(string(observer,format='(A)'),2) endif else obsName='500' ;--I.7-- Retrieve KOALA Tools Variables ------------------------------------------------------ confAstroIM= initIDL() confKOALA = initIDL(confAstroIM.soft.koala, /KOALA) ;--I.X-- Useful Constants -------------------------------------------------------------------- c_kms = 299792.458D ;-km/s AUtoKM = 149598000.D ;-km c_aus = c_kms / AUtoKM tempoIN = '/tmp/lc4generation.dat' tempoLC = '/tmp/lc.dat' ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- II -- Generate Lightcurve from 3-D Shape -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--II.1-- Ephemerides of the Sun ------------------------------------------------------------- sunEph = voMiriade_callEphemCC(id,ep,/Helio,/Julian,tcoor=2,rPlane=2,nbd=nbd,step=step,/Web,dump='/tmp/sun') sunEph = sunEph[0:n_elements(sunEph)-2] ;--II.2-- Ephemerides of the Observer -------------------------------------------------------- obsEph = voMiriade_callEphemCC( id, ep, observer=obsName, /Julian, tcoor=2, rPlane=2, $ nbd=nbd, step=step, /Web,dump='/tmp/obs' ) obsEph = obsEph[0:n_elements(obsEph)-2] ;--II.3-- Light-Time Correction -------------------------------------------------------------- LTC = (sqrt(obsEph.dObs) / c_aus)/3600./24.D jdLTC = obsEph.JD - LTC ;--II.4-- Lightcurve Generation -------------------------------------------------------------- forprint, jdLTC, fltarr(nbd), $ -sunEph.x, -sunEph.y, -sunEph.z, $ -obsEph.x, -obsEph.y, -obsEph.z, $ comment=transpose(['1',strtrim(string(nbd),2)+' 0']), $ format='(D16.8,2x,F4.1,2x,6(2x,F11.6))', textout=tempoIn, /Silent spawn, 'cat '+tempoIn+$ ' | '+confKOALA.lc.concav+' -i '+$ spin+' '+$ confKOALA.lc.config+' '+$ shape+' > '+tempoLC, genO, genE lc = koalaReadLC( tempoLC ) *lc.jd = obsEph.jd ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- III -- Generate Period-Long Lightcurve for Normalization -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--III.1-- Avoid Computation if LC requested for Period -------------------------------------- if keyword_set(period) then begin meanFlux = mean( *lc.flux ) ;--III.2-- Generate the LC over the Period --------------------------------------------------- endif else begin ;--III.2.1-- Set nbd & step ---------------------------------------------------------------- nbdNorm = 720 stepNorm = strtrim(string(spinStr.ph/nbdNorm,format='(F11.6)'),2)+'h' ;--III.2.2-- Ephemerides of the Sun -------------------------------------------------------- sunNorm = voMiriade_callEphemCC(id,ep,/Helio,/Julian,tcoor=2,rPlane=2,nbd=nbdNorm,step=stepNorm,/Web) sunNorm = sunNorm[0:n_elements(sunNorm)-2] ;--III.2.3-- Ephemerides of the Observer --------------------------------------------------- obsNorm = voMiriade_callEphemCC( id, ep, observer=obsName, /Julian, tcoor=2, rPlane=2, $ nbd=nbdNorm, step=stepNorm, /Web ) obsNorm = obsNorm[0:n_elements(obsNorm)-2] ;--III.2.4-- Light-Time Correction --------------------------------------------------------- normLTC = (sqrt(obsNorm.dObs) / c_aus)/3600./24.D jdNorm = sunNorm.JD - normLTC ;--III.2.5-- Lightcurve Generation --------------------------------------------------------- forprint, jdNorm, fltarr(nbdNorm), $ -sunNorm.x, -sunNorm.y, -sunNorm.z, $ -obsNorm.x, -obsNorm.y, -obsNorm.z, $ comment=transpose(['1',strtrim(string(nbdNorm),2)+' 0']), $ format='(D16.8,2x,F4.1,2x,6(2x,F11.6))', textout=tempoIn, /Silent spawn, 'cat '+tempoIn+$ ' | '+confKOALA.lc.concav+' -i '+$ spin+' '+$ confKOALA.lc.config+' '+$ shape+' > '+tempoLC, genO, genE lcNorm = koalaReadLC( tempoLC ) meanFlux = mean(*lcNorm.flux) normMag = 2.5*alog10( *lcNorm.flux ) endelse ;--III.3-- Normalize the Flux of the Lightcurve ---------------------------------------------- *lc.flux /= meanFlux ; print, 'meanflux', meanFlux ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- IV -- Convert Flux from 3-D shape into Apparent Magnitude ---------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--IV.1-- Convert Flux to Magnitude ---------------------------------------------------------- lcMag = -2.5*alog10(*lc.flux) ;--IV.2-- Compute Phase Function ------------------------------------------------------------- ;--IV.2.1-- Sun-Target and Obs-Target Vectors vSun = [[sunEph.X], [sunEph.Y], [sunEph.Z]] vObs = [[obsEph.X], [obsEph.Y], [obsEph.Z]] ;--IV.2.2-- Sun-Target and Obs-Target Distances dSun = sqrt( total(vSun*vSun,2) ) dObs = sqrt( total(vObs*vObs,2) ) ;--IV.2.3-- Phase Angle phase = acos( total(vSun*vObs,2)/(dSun*dObs) ) ;--IV.2.4-- Bowell Phase Function G = 0.15 phi1 = exp( -3.33 *tan(phase/2.)^0.63 ) phi2 = exp( -1.87 *tan(phase/2.)^1.22 ) phaseFunction = -2.5*alog10( (1-G)*phi1 + G*phi2 ) ;--IV.3-- Compute Distance Function ---------------------------------------------------------- distFunction = 2.5*alog10( dObs*dObs * dSun*dSun ) ;--IV.4-- Retrieve Target's Absolute Magnitude ----------------------------------------------- ssoInfo = astElements( id ) ;--IV.5-- Compute Target's Apparent Magnitude ------------------------------------------------ mag = ssoInfo.H + distFunction + phaseFunction + lcMag ; sep = fltarr( nbd, 2) ; ssp = fltarr( nbd, 2) ; for k=0, nbd-1 do begin ; eco = subObserver_coord( -vObs[k,*], spinStr, (*lc.jd)[k]) ; sco = subObserver_coord( -vSun[k,*], spinStr, (*lc.jd)[k]) ; sep[k,0] = eco.lon ; sep[k,1] = eco.lat ; ssp[k,0] = sco.lon ; ssp[k,1] = sco.lat ; endfor ; ; eph=voMiriade_callEphemCC(id, ep,nbd=nbd,step=step,/web) ; eph = eph[0:n_elements(eph)-2] ; print, 'h', ssoInfo.h ; print, 'dist', median( distFunction ) , max(distFunction )-min(distFunction ) ; print, 'phase', median( phaseFunction ), max(phaseFunction)-min(phaseFunction) ; print, 'cdr', median(lcMag), max( lcMag ) - min(lcMag) ; print, 'mag', median(mag), max( Mag ) - min(Mag) ; print, 'eph', median(eph.mag), max( eph.Mag ) - min(eph.Mag) ; ; ; cgPlot, /yNoZero, *lc.jd, mag, color='Cornflower blue',yr=[10,18], $ ; xtitle='Time (day)', yTitle='Apparent Magnitude' ; cgPlot, /OverPlot, *lcNorm.jd, ssoInfo.H + distFunction + phaseFunction + normMag, color='Orange' ; cgPlot, /OverPlot, eph.jd, eph.mag, color='Red' ; cgPlot, /yNoZero, sunNorm.JD, normMag, color='Cornflower blue',yr=[10,12], $ ; xtitle='Time (day)', yTitle='Apparent Magnitude' ; cgPlot, /OverPlot, *lc.jd, lcMag, color='Red', psym='Filled Circle' ;stop ; cgPlot, /OverPlot, *lcNorm.jd, ssoInfo.H + distFunction + phaseFunction + normMag, color='Orange' ; cgPlot, /OverPlot, eph.jd, eph.mag, color='Red' ; ; ; cgPlot, eph.jd, eph.mag - (ssoInfo.H+phasefunction+distfunction), $ ; xtitle='Time (day)', yTitle='Miriade - (mon Bowell)' ; ;stop ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- V -- Export & Disk I/O ---------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; case format of ;--V.1-- Simple: JD vs Mag ------------------------------------------------------------------- 1: begin if keyword_set(dump) then $ forprint, *lc.JD, mag, comment='JulianDay, Mag', $ format='(D16.8,",",F5.2)', textout=dump, /Silent out = replicate({jd:0.d,mag:0.}, nbd) out.jd = *lc.jd out.mag = mag return, out end ;--V.2-- Simple: JD vs dMag ------------------------------------------------------------------ 2: begin if keyword_set(dump) then $ forprint, *lc.JD, lcMag, comment='JulianDay, Mag', $ format='(D16.8,",",F5.2)', textout=dump, /Silent out = replicate({jd:0.d,mag:0.}, nbd) out.jd = *lc.jd out.mag = lcMag return, out end ;--V.3-- Simple: JD vs Flux ------------------------------------------------------------------ ;--V.X-- Else... ----------------------------------------------------------------------------- else: begin print, 'code the test at the beginning!!' end endcase ; ;--III.1-- Write the Lightcurve Input File on Disk ; if keyword_set(dump) then $ ; forprint, JD, fltarr(nbd), $ ; -sunEph.x, -sunEph.y, -sunEph.z, $ ; -obsEph.x, -obsEph.y, -obsEph.z, $ ; comment=transpose(['1',strtrim(string(nbd),2)+' 0']), $ ; format='(D16.8,2x,F4.1,2x,6(2x,F11.6))', $ ; textout=dump, /SILENT ; ; ;-III.2-- Concatenate Information into a Structure ; RETURN, {num: nbd, jd: JD, $ ; sun: -[ [sunEph.x], [sunEph.y], [sunEph.z] ], $ ; obs: -[ [obsEph.x], [obsEph.y], [obsEph.z] ] } end