; ; ; simulate time-resolved photometry of euclid VIS/NISP ; ; use taxonomy colors ; randomize LC parameters = [amplitude, period] ; ; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- I -- Initialization And Input Verification -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; doFigOnly = 01 if doFigOnly eq 1 and keyword_set(uncP_vs_A) then goto, j2Fig ;--I.1-- Define Directories ------------------------------------------------------------------; spawn, 'hostname', host case strTrim(host) of 'hyperion': dirEuclid = '/data/euclid/' 'endymion': dirEuclid = '/home/bcarry/work/data/euclid/' else: stop endcase dirLC = dirEuclid + 'lc/' dirPHOT = dirLC+'simuPHOT/' dirTaxo = dirEuclid + 'taxo/' ;--I.2-- Simulation Parameters ---------------------------------------------------------------; ;--I.2.1-- Period and Amplitude Range ; period={ min:1., max:200., N:50L } ; amp ={ min:0., max:1.6, N:8L } period={ min:1., max:200., N:100l} amp ={ min:0., max:1.6, N:8L } nbLC = amp.N * period.N ;--I.2.2-- Period and Amplitude Histogram binA = 0.1 nbBinA = ( amp.max - amp.min) / binA ; binP = 0.25 ; HOUR nbBinP = period.N ;(period.max - period.min) / binP vecP = cgLogGen( period.N, start=period.min, finish=period.max ) ;--I.2.3-- Noise seed=1 sigNoise = 0.02 ;--I.2.4-- Residuals and Uncertainty res = { min:-0.2, max:0.2, bin:0.005, N:0L } unc = { min:0., max:0.2, bin:0.005, N:0L } res.N = (res.max-res.min)/res.bin unc.N = (unc.max-unc.min)/unc.bin ;--I.3-- Read Bus-DeMeo Asteroid Colors ------------------------------------------------------; readcol, dirTaxo+'euclid-ssomag.csv', delimiter=',', /Silent, $ ssoNum, ssoName, ssoClass, $ VmY, VmJ, VmH, YmJ, YmH, JmH, $ dVmY, dVmJ, dVmH, dYmJ, dYmH, dJmH, $ classId, compLabel, complexId, $ format='(L,A,A,'+$ 'F,F,F,F,F,F,'+$ 'F,F,F,F,F,F,'+$ 'I,A,I)' ; VmY = 0.9 ; VmJ = 1.1 ; VmH = 1.4 ; YmJ = VmJ - VmY ; YmH = VmH - VmY ; JmH = VmH - VmJ ; ssoNum=0 ; ssoName='toto' nbSSO = n_elements(VmY) ; nbSSO = 1 ssoName=strTrim(ssoName,2) ;--I.4-- Read Merlin KBO Colors --------------------------------------------------------------; ;TBD: ; 1) check euclid-vnir to export KBO colors 1-by-1 ; 2) read it like BDM ; ;--I.5-- Time Variables ----------------------------------------------------------------------; time = findgen(4500) ; 1h+ long LC in SECONDS ;-- and convert period in seconds!!! period.min *= 3600. period.max *= 3600. vecP *= 3600 ;--I.6-- Sequence of observation -------------------------------------------------------------; ;--I.6.1-- Read Sequence of Observation readcol, dirPHOT+'time.def', delimiter=',', /Silent, $ filtId, obsStart, obsDIT, format='(A,F,F)' nbObs = n_elements(filtId) filtId=strTrim(filtId,2) ;--I.6.2-- Check for Inconsistency outtaRange = where( time le min(obsStart) or time ge max(obsStart+obsDIT), nbRange ) if nbRange lt 2 then begin print, 'Observing sequence outside time variable: ' print, ' Time: '+strTrim(string(min(time)),2)+' -- '+ strTrim(string(max(time)),2) print, ' Obs: ' +strTrim(string(min(obsStart)),2)+' -- '+strTrim(string(max(obsStart+obsDIT)),2) stop endif ;--I.6.3-- Identify Set of Filters -------------------------------------------------------; curFilt = filtId uniqFilt = uniq(curFilt, sort(curFilt)) listFilt = curFilt[uniqFilt] listFilt = ['VIS','Y','J','H'] nbUniqFilt = n_elements( uniqFilt ) colFilt=['Black', 'Green', 'Cornflower Blue', 'Red'] symFilt=['Filled Circle', 'Open Square', 'Open Circle', 'OpenUpTriangle'] nbComb = 0 for kF=0, nbUniqFilt-1 do for kO=kF+1, nbUniqFilt-1 do nbComb++ ;--I.7-- Giantic Color Variable --------------------------------------------------------------; nbMeth = 2 stat = replicate( {period:0., amp:0., $ ;-LC parameters col:fltarr(nbComb, nbMeth), $ ;-Average colors unc:fltarr(nbComb, nbMeth), $ ;-Standard deviation colors res:fltarr(nbComb, nbMeth), $ ;-Residuals vs Simulated N:fltarr(nbComb, nbMeth)}, $ ;-Number of measures nbSSO, nbLC ) ;--I.8-- Distribution of SSOs in Period vs Amplitude -----------------------------------------; lcPDS = pds_lcparam_getcatalog(/dump) val = where( strCmp(lcPDS.p.qual,'2',1) or strCmp(lcPDS.p.qual,'3',1) ) lcPDS = lcPDS[val] ; pds2D = hist_2d( lcPDS.a.max, lcPDS.p.val*3600., $ ; min1=amp.min, max1=amp.max, bin1=binA, $ ; min2=period.min, max2=period.max, bin2=binP ) ; pdsPDF = pds2D / total(pds2d) t0 = systime(/JULIAN,/UTC) bigLC = fltarr(nbSSO, 4500, nbLC ) mag = fltarr(nbSSO, nbLC, nbObs) ;--I.9-- Loop Over Input SSO Taxonomies ------------------------------------------------------; for kSSO=0, nbSSO-1 do begin ;--I.9.1-- Build Reference Filter-Filter Matrix for Current SSO refCol = fltarr(nbUniqFilt,nbUniqFilt) refCol[0,*] = [ 0 , VmY[kSSO], VmJ[kSSO], VmH[kSSO] ] refCol[1,*] = [-VmY[kSSO], 0 , YmJ[kSSO], YmH[kSSO] ] refCol[2,*] = [-VmJ[kSSO], -YmJ[kSSO], 0 , JmH[kSSO] ] refCol[3,*] = [-VmH[kSSO], -YmH[kSSO], -JmH[kSSO], 0 ] ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- II -- Compute Synthetic Photometry -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--II.1-- Generate Random Lightcurves ------------------------------------------------------; ;--II.1.1-- Blind Generation lcs = randomLC(time, seed+kSSO, pMin=period.min, pMax=period.max, nbP=period.N, /pLog, $ aMin=amp.min, aMax=amp.max, nbA=amp.N, oPeriod=oPeriod, oAmplitude=oAmplitude ) bigLC[kSSO,*,*]=lcs ;--II.1.2-- Store LC Parameters: Period and Amplitude for kLC=0, nbLC-1 do begin stat[kSSO,kLC].period = oPeriod[kLC] stat[kSSO,kLC].amp = oAmplitude[kLC] endfor ;--II.2-- Simulate Photometry --------------------------------------------------------------; noise = randomn(seed+4, nbLC, nbObs) * sigNoise for kO=0, nbObs-1 do begin ;--II.2.1-- Time Range of Current Filter range = where( time ge obsStart[kO] and time le (obsStart[kO]+obsDIT[kO])) ;--II.2.2-- Current Filter Color case filtId[kO] of 'VIS': locColor = 0. 'Y': locColor = VmY[kSSO] 'J': locColor = VmJ[kSSO] 'H': locColor = VmH[kSSO] endcase ;--- tbd store reference Fi-Fj colors in stat structure ;--II.2.3-- Photometry *WITH SED COLORS* for kLC=0, nbLC-1 do begin mag[kSSO, kLC,kO] = mean(lcs[range,kLC]) - locColor + noise[kLC,kO] ; cgPlot, time, lcs[*,kLC] ; cgPlot, /OverPlot, time[range], lcs[range,kLC], color='red', thick=4 ; cgPlot, /OverPlot, median(time[range]), median(lcs[range,kLC]), color='blue', psym='Filled Circle' ; stop endfor endfor ;--II.2-- End of loop on Observing Time ;--II.4-- Write Photometry to Disk ---------------------------------------------------------; nameLC=string(ssoNum[kSSO],format='(I06)')+'-'+ssoName[kSSO]+'.csv' openW, idEx, dirPHOT+'lc/'+nameLC, /Get_Lun printf, idEx, 'lcID, Time, Mag, Filt' for kLC=0, nbLC-1 do for kF=0, nbObs-1 do begin printf, idEX, kLC, (obsStart[kF]+obsDIT[kF]/2.), mag[kSSO, kLC,kF], filtId[kF], $ format='(I6,",",F7.2,",",F9.5,",",A-3)' endfor free_lun, idEx ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- III -- Determine Colors from Photometry -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; for kLC=0, nbLC-1 do begin ;--III.1-- Set Up Local Variables --------------------------------------------------------; curLC = lcs[*,kLC] curTime = obsStart+obsDIT/2. curMag = mag[kSSO, kLC,*] curFilt = filtId ; window, 2 ; cgPlot, time, curLC, yRange=[-1,1]*amp.max, xTitle='Time (s)', yTitle='Mag' ; cgPlot, /OverPlot, time, curLC-VmY[kSSO], color='Red' ; cgPlot, /OverPlot, time, curLC-VmJ[kSSO], color='Green' ; cgPlot, /OverPlot, time, curLC-VmH[kSSO], color='Navy' ; cgPlot, /OverPlot, curTime, curMag, psym='Filled Circle', color='Blue' colV = fltarr( nbMeth, nbUniqFilt, nbUniqFilt, nbObs ) colN = intarr( nbMeth, nbUniqFilt, nbUniqFilt ) ;--III.3-- Trial of Methods to Measure Colors --------------------------------------------; ;--III.3.1-- Method 0: Interpolation over N Points ---------------------------------------; kMeth=0 for kF=0, nbUniqFilt-1 do begin ;--III.3.1/A-- Main Filter Id fId = where( strCmp(curFilt, listFilt[kF]), nbLocFilt ) for kLF=0, nbLocFilt-2 do begin locArr = fID[ kLF:kLF+1 ] ; cgPlot, /OverPlot, curTime[locArr], curMag[locArr], psym='Open Circle', color='Red', symSize=2+kLF ;--III.3.1/B-- Loop Over Observations In-Between bracket = where( curTime gt curTime[locArr[0]] and curTime lt curTime[locArr[1]], nbBracket ) for kB=0, nbBracket-1 do begin ;--III.3.1/C-- Find other Filter Id oFilt = curFilt[bracket[kB]] oId = where( strCmp(oFilt, listFilt) ) oId=oId[0] ;--III.3.1/D-- Find other Filter Mag + Time oTime = curTime[bracket[kB]] oMag = curMag[bracket[kB]] ;--III.3.1/E-- Interpolate interp = interpol( curMag[fId], curTime[fId], oTime, /Quadratic ) colV[kMeth, kF, oId, colN[kMeth,kF,oId]] = interp-oMag colN[kMeth, kF, oId]++ ; cgPlot, /OverPlot, [oTime,oTime], [oMag,interp] ; cgPlot, /OverPlot, oTime, oMag, psym='Open Square', color='Green', symSize=2+kLF ; cgPlot, /OverPlot, oTime, interp, psym='Open circle', color='black' ; cgText, oTime+50, oMag, oFilt, color='Green' endfor ;--III.3.1/B-- End of Loop on Bracket endfor ;--End of Loop on Filter Times endfor ;--III.3.1-- End of Loop on Filters ;--III.3.2-- Method 1: Color with Closest in Time ----------------------------------------; kMeth=1 for kF=0, nbUniqFilt-1 do begin ;--III.3.2/A-- Main Filter Id fId = where( strCmp(curFilt, listFilt[kF]), nbLocFilt ) for kLF=0, nbLocFilt-1 do begin locArr = fID[kLF] ; cgPlot, /OverPlot, curTime[locArr], curMag[locArr], psym='Open Circle', color='Red', symSize=2+kLF ;--III.3.2/B-- Loop Over Other Filters for kO=0, nbUniqFilt-1 do if kO ne kF then begin oId = where( strCmp(curFilt, listFilt[kO]) ) ; cgPlot, /OverPlot, curTime[oId], curMag[oId], psym='Open Square', color='Blue', symSize=2+kLF ;--III.3.2/C-- Find Closest in Time deltaT = abs( curTime[oId] - curTime[locArr]) minT = min(deltaT, pMinDeltaT ) oTime = curTime[oId[pMinDeltaT]] oMag = curMag[ oId[pMinDeltaT]] ;--III.3.2/D-- Simple Magnitude Difference colV[kMeth, kF, kO, colN[kMeth,kF,kO]] = curMag[locArr]-oMag colN[kMeth, kF, kO]++ ; cgPlot, /OverPlot, oTime, oMag, psym='Open Square', color='Orange', symSize=2+kLF endif endfor ;--End of Loop on Filter Times endfor ;--III.3.2-- End of Loop on Filters ;--III.4-- Compile All Filter-Filter Measurements ----------------------------------------; for kMeth=0, nbMeth-1 do begin kComb=0 for kF=0, nbUniqFilt-1 do for kO=kF+1, nbUniqFilt-1 do begin ;--III.4.1-- Filter1 - Filter2 nbEval1 = colN[kMeth,kF,kO] if nbEval1 gt 0 then evalArr1 = colV[kMeth, kF,kO,0:nbEval1-1] ;--III.4.2-- Filter2 - Filter1 nbEval2 = colN[kMeth,kO,kF] if nbEval2 gt 0 then evalArr2 = - colV[kMeth, kO,kF,0:nbEval2-1] ;--III.4.3-- Number of Evaluations nbEval = nbEval1 + nbEval2 if nbEval gt 0 then begin if nbEval1 gt 0 and nbEval2 gt 0 then evalArr=[evalArr1,evalArr2] else begin if nbEval1 gt 0 then evalArr = evalArr1 if nbEval2 gt 0 then evalArr = evalArr2 endelse ;--III.4.4-- Number, Mean, Stddev of Color: F1-F2 stat[kSSO,kLC].N[kComb,kMeth] = nbEval stat[kSSO,kLC].col[kComb,kMeth] = mean(evalArr) stat[kSSO,kLC].res[kComb,kMeth] = stat[kSSO,kLC].col[kComb,kMeth] - refCol[kF,kO] stat[kSSO,kLC].unc[kComb,kMeth] = stddev(evalArr) kComb++ endif endfor ;--III.4-- End of Loop on Filter Combinations endfor ;--III.4-- End of Loop on Methods endfor ;--II-- End of Loop on LC ;stop endfor ;--I.8-- End of Loop on SSOs t1 = systime(/JULIAN,/UTC) print, 'Simulation of photometry finished: '+string( (t1-t0)*24*60, format='(F5.2)') ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- IV -- Compute Average Colors & Residuals -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--IV.1-- Average Color Measurements ---------------------------------------------------------; openW, idEx, dirPHOT+'lc/averagecolor-'+string(sigNoise,format='(F5.3)')+'.csv', /Get_Lun printf, idEx, 'methId, colorId, N, color, dColor, SNR, residual, dResidual' for kM=0, nbMeth-1 do begin print, '' print, ' Color N Color Uncert SNR meanRes stdRES' kComb=0 for kF=0, nbUniqFilt-1 do for kO=kF+1, nbUniqFilt-1 do begin colId = listFilt[kF]+'-'+listFilt[kO] meanCol = meanWithUnc( stat.col[kComb,kM], stat.unc[kComb,kM] ) meanRes = mean( stat.res[kComb,kM] ) stdRes = stddev( stat.res[kComb,kM] ) nbEval = total( stat.N[kComb,kM] ) kComb++ print, kM, colId, nbEval, meanCol[0], meanCol[1], meanCol[1]/sigNoise, $ meanRes, stdRes, $ format='(I1,1x,A-8,2x,I7,2x,F7.4,2x,F7.4,2x,F7.3,2x,F7.4,2x,F7.4)' printf, idEx, kM, colId, nbEval, meanCol[0], meanCol[1], meanCol[1]/sigNoise, meanRes, stdRes, $ format='(I1,", ",A-8,", ",I7,", ",F7.4,", ",F7.4,", ",F7.3,", ",F7.4,", ",F7.4)' endfor endfor free_lun, idEx ;--IV.2-- Build Map of Residuals vs Parameters -----------------------------------------------; colP_vs_A = fltarr(nbBinA, nbBinP, nbMeth, nbComb) resP_vs_A = fltarr(nbBinA, nbBinP, nbMeth, nbComb) uncP_vs_A = fltarr(nbBinA, nbBinP, nbMeth, nbComb) pds2D = fltarr(nbBinA, nbBinP ) for kA=0, nbBinA-1 do begin for kP=0, nbBinP-2 do begin ;-Euclid Data cur = where( stat.period ge vecP[kP] and stat.period lt vecP[kP+1] and $ stat.amp ge kA*binA and stat.amp lt (kA+1)*binA, nbCur ) if cur[0] ne -1 then begin for kM=0, nbMeth-1 do begin for kC=0, nbComb-1 do begin colP_vs_A[kA,kP,kM,kC] = mean( stat[cur].col[kC,kM], /NaN ) resP_vs_A[kA,kP,kM,kC] = mean( stat[cur].res[kC,kM], /NaN ) uncP_vs_A[kA,kP,kM,kC] = mean( stat[cur].unc[kC,kM], /NaN ) endfor endfor endif ;-PDS data cur = where( lcPDS.p.val*3600. ge vecP[kP] and lcPDS.p.val*3600. lt vecP[kP+1] and $ lcPDS.a.max ge kA*binA and lcPDS.a.max lt (kA+1)*binA, nbCur ) if cur[0] ne -1 then pds2D[kA,kP] = nbCur endfor ;--IV.2-- End of Loop on Color Combinaisons endfor ;--IV.2-- End of Loop on Methods pdsPDF = pds2d / total(pds2d) t2 = systime(/JULIAN,/UTC) print, 'Compilation of results finished: '+string( (t2-t1)*24*60, format='(F5.2)') ;--IV.3-- PDF of Residuals or Uncertainty ----------------------------------------------------; ;- resArr = findgen(res.N)*res.bin ;- resPDF = fltarr(res.N, nbMeth, nbComb) ;- kM=0 ;- kC=0 ;- for kR=0, res.N-2 do begin ;- ;- loc = resP_vs_A[*,*,kM,kC] ;- sel = where( loc ge resArr[kR] and loc lt resArr[kR+1] ) ;- if loc[0] ne -1 then begin ;- resPDF[kR,kM,kC] += total(pdsPDF[sel]) ;- endif ;- endfor print, '' print, 'Residuals' print, 'mean 0, mean1, std1, std2, ratioMean, ratioStd' for kC=0, nbComb-1 do $ print, mean( resP_vs_A[*,*,0,kC]) , mean( resP_vs_A[*,*,1,kC]), $ stddev( resP_vs_A[*,*,0,kC]) ,stddev( resP_vs_A[*,*,1,kC]), $ mean( resP_vs_A[*,*,0,kC]) /mean( resP_vs_A[*,*,1,kC]), $ stddev( resP_vs_A[*,*,1,kC])/stddev( resP_vs_A[*,*,0,kC]) print, '' print, 'Dispersion' print, 'mean 0, mean1, std1, std2, ratioMean, ratioStd' for kC=0, nbComb-1 do $ print, mean( uncP_vs_A[*,*,0,kC]) , mean( uncP_vs_A[*,*,1,kC]), $ stddev( uncP_vs_A[*,*,0,kC]) ,stddev( uncP_vs_A[*,*,1,kC]), $ mean( uncP_vs_A[*,*,0,kC]) /mean( uncP_vs_A[*,*,1,kC]), $ stddev( uncP_vs_A[*,*,1,kC])/stddev( uncP_vs_A[*,*,0,kC]) ; cghistoplot, resP_vs_A[*,*,1,0] ; cghistoplot, resP_vs_A[*,*,0,0],/oplot ; ; cghistoplot, uncP_vs_A[*,*,1,0] ; cghistoplot, uncP_vs_A[*,*,0,0],/oplot writefits, dirPHOT+'resP_vs_A.fits', resP_vs_A writefits, dirPHOT+'uncP_vs_A.fits', uncP_vs_A writefits, dirPHOT+'colP_vs_A.fits', colP_vs_A openW, idCol, dirPHOT+'colP_vs_A.csv', /Get_Lun printf, idCol, 'Amplitude, Period, Method, Color, Value' openW, idRes, dirPHOT+'resP_vs_A.csv', /Get_Lun printf, idRes, 'Amplitude, Period, Method, Color, Residual' openW, idUnc, dirPHOT+'uncP_vs_A.csv', /Get_Lun printf, idUnc, 'Amplitude, Period, Method, Color, Dispersion' print, nbMeth, nbComb, nbBinA, nbBinP help, resP_vs_A vecPout=vecP for kP=0, nbBinP-2 do vecPout[kP] = mean(vecP[kp:kP+1])/3600. kM=0 & kC=0 & kA=0 & kP=0 for kM=0, nbMeth-1 do begin for kC=0, nbComb-1 do begin for kA=0, nbBinA-1 do begin for kP=0, nbBinP-2 do begin printf, idRes, (kA+0.5)*binA, vecPout[kP], kM, kC, resP_vs_A[kA,kP,kM,kC], $ format='(F4.2,",",F8.2,",",I1,",",I1,",",F9.4)' printf, idUnc, (kA+0.5)*binA, vecPout[kP], kM, kC, uncP_vs_A[kA,kP,kM,kC], $ format='(F4.2,",",F8.2,",",I1,",",I1,",",F9.4)' printf, idCol, (kA+0.5)*binA, vecPout[kP], kM, kC, colP_vs_A[kA,kP,kM,kC], $ format='(F4.2,",",F8.2,",",I1,",",I1,",",F9.4)' endfor endfor endfor endfor free_lun, idCol free_lun, idRes free_lun, idUnc ; print, mean( resP_vs_A[*,*,0,kC])/stddev( resP_vs_A[*,*,0,kC]), $ ; mean( resP_vs_A[*,*,1,kC])/stddev( resP_vs_A[*,*,1,kC]) ; t3 = systime(/JULIAN,/UTC) print, 'Simulations finished: '+string( (t3-t0)*24*60, format='(F5.2)') ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- V -- Figures Histograms -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; j2Fig: set_plot, 'X' ;--V.1-- Compute Histograms ------------------------------------------------------------------; binSig = 0.001 binRes = 0.001 case sigNoise of ;--- Good 0.020 : begin binSig = 0.0005 sigRang = [0.016,0.024] resRang = [-0.01,0.01] xTickInt1=0.002 xTickInt2=0.005 end ;--Median 0.05 : begin sigRang = [0.040,0.100] resRang = [-0.03,0.03] xTickInt1=0.02 xTickInt2=0.010 end ;--- Poor 0.10 : begin sigRang = [0.08,0.15] resRang = [-0.03,0.03] xTickInt1=0.02 xTickInt2=0.010 end endcase nbTickSig = (sigRang[1] - sigRang[0])/binSig + 1 nbTickRes = (resRang[1] - resRang[0])/binRes + 1 cgHistoPlot, uncP_vs_A[*,*,0,0], loc=xSig, histData=ySig, binSize=binSig, minInput=sigRang[0], maxInput=sigRang[1] nbSig = n_elements(ySig) hSig = fltarr( nbSig, 2, 3 ) cgHistoPlot, resP_vs_A[*,*,0,0], loc=xRes, histData=yRes, binSize=binRes, minInput=resRang[0], maxInput=resRang[1] nbRes = n_elements(yRes) hRes = fltarr( nbRes, 2, 3 ) for kM=0, 1 do begin for kC=0, 3-1 do begin cgHistoPlot, uncP_vs_A[*,*,kM,kC], loc=xSig, histData=ySig, binSize=binSig, minInput=sigRang[0], maxInput=sigRang[1] hSig[*,kM,kC] = ySig cgHistoPlot, resP_vs_A[*,*,kM,kC], loc=xRes, histData=yRes, binSize=binRes, minInput=resRang[0], maxInput=resRang[1] hRes[*,kM,kC] = yRes endfor endfor wdelete ;--V.2-- EPS/PNG Output ----------------------------------------------------------------------; cgPS_open, xSize=22, ySize=26.5, Filename=dirPHOT+'euclid-acolr-'+string(sigNoise,format='(F5.3)')+'.eps', $ /Metric, /Decomposed, /Landscape, /Encapsulated, Language_Level=2, /Quiet xBot = 2200. xLen = 9000. yBot = 2000. yLen = 8000. maxH = 1.1 yRang=[0,maxH] colCum='Cornflower Blue' for kC=0, 3-1 do begin ;--V.3-- Axes 'n' Everything -----------------------------------------------------------------; posArr1 = [xBot, yBot+kC*yLen, xBot+xLen, yBot+(kC+1)*yLen] posArr2 = [xBot+xLen, yBot+kC*yLen, xBot+xLen*2., yBot+(kC+1)*yLen] if kC eq 0 then begin xTitle1='Dispersion (mmag)' xTitle2='Residuals (mmag)' xTickN1= [strTrim(string( (sigRang[0] + xTickInt1*findgen(nbTickSig-1))*1000, format='(I3)'),2),' '] xTickN2= [strTrim(string( (resRang[0] + xTickInt2*findgen(nbTickRes))*1000, format='(I3)'),2)] endif else begin xTitle1='' xTitle2='' xTickN1= replicate(' ',15) xTickN2= replicate(' ',15) endelse ;--V.4-- Dispersion of Measurements ----------------------------------------------------------; cgPlot, sigRang, [0,maxH], /NoData, /NoErase, /Device, $ xStyle=1, xRange=sigRang, xTickInt=xTickInt1, xMinor=5, $ yStyle=1, yRange=yRang, yTickInt=0.2, yMinor=5, $ xTickLen=0.04, $ yTickLen=0.03, $ charSize=1.5, position=posArr1, $ xTitle=xTitle1, xTickName=xTickN1, $ yTitle='Distribution ('+listFilt[0]+'-'+listFilt[kC+1]+')' cgPlot, /OverPlot, xSig, hSig[*,0, kC]/max(hSig[*,0,kC]), psym=10 cgPlot, /OverPlot, xSig, total(hSig[*,0, kC],/Cum)/total(hSig[*,0,kC]), color=colCum ;--V.5-- Residuals of Measurements -----------------------------------------------------------; cgPlot, resRang, [0,maxH], /NoData, /NoErase, /Device, $ xStyle=1, xRange=resRang, xTickInt=xTickInt2, xMinor=5, $ yStyle=9, yRange=yRang, yTickInt=0.2, yMinor=5, $ xTickLen=0.04, $ yTickLen=0.03, $ charSize=1.5, position=posArr2, $ xTitle=xTitle2, xTickName=xTickN2, $ yTitle='', yTickName=replicate(' ',20) cgPlot, /OverPlot, xRes, hRes[*,0 ,kC]/max(hRes[*,0,kC]), psym=10 cgPlot, /OverPlot, xRes, total(hRes[*,0 ,kC],/Cum)/total(hRes[*,0,kC]), color=colCum cgAxis, /YAxis, yStyle=1, yRange=yRang, color=colCum, resRang[1], charSize=1.5, yTitle='Cumulative distribution' endfor cgPS_close, /png, Delete_PS=0 ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- V -- Figures 2d Distro of Results -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; cgPS_open, xSize=22, ySize=26.5, Filename=dirPHOT+'euclid-p_vs_a-'+string(sigNoise,format='(F5.3)')+'-res.eps', $ /Metric, /Decomposed, /Landscape, /Encapsulated, Language_Level=2, /Quiet ;--V.1-- Parameters --------------------------------------------------------------------------; xBot = 0.15 xMid = 0.8 xTop = 0.9 yBot = 0.1 yTop = 0.95 pArr1 = [ xBot, yBot, xMid, yTop ] pArr2 = [ xMid, yBot, xTop, yTop ] cgLoadCT, 39, NCOLORS=255, RGB_TABLE=palette ;--V.3-- Residuals ---------------------------------------------------------------------------; resImg = cgImgScl( resP_vs_A[*,*,0,0], minValue=resRang[0], maxValue=resRang[1] ) cgImage, /NoErase, resImg, /Normal, position=pArr1, palette=palette lvlArr = [0.00015, 0.00283, 0.0056] contour, pdsPDF, min_value=min(pdsPDF), max_value=max(pdsPDF), $ levels=lvlArr, Path_Info=path, Path_Xy=XY, /Closed, /Isotropic, /Path_Data_Coords nbCurve=n_elements(path) colSSO = 'White' yShi=[10, 3, 0.5] for kP=0,1 do begin pathIndex= [indgen(path(kP).N), 0] c = xy(*,path(kP).offset + pathIndex ) cgPlot, /OverPlot, c[0,*], c[1,*], color=colSSO ww= where( pdsPDF ge lvlArr[kP]) frac = string(total(pdsPDF[ww])*100.,format='(I3)') print, 'PDS : '+string(lvlArr[kP],format='(F5.3)')+' : '+frac+' %' yL = max(c[1,*],pMax) xL = c[0,pMax] cgText, /Data, xL, yL, frac+'%', color=colSSO endfor cgPlot, 0,vecP/3600., /NoData, /NoErase, /Normal, $ position=pArr1, $ xTickLen=0.02, yTickLen=0.04, $ xStyle=1, xRange=[amp.min, amp.max], $ yStyle=1, yRange=[period.min, period.max]/3600., /yLog, $ xTitle='Lightcurve amplitude (mag)', $ yTitle='Rotation period (h)' xKey=[0.8,1.1] yKey=130*[1,1] cgPlot, xKey, yKey, /OverPlot, color=colSSO cgText, /Data, 1.3, yKey+5, 'Known' , align=0.5, color=colSSO cgText, /Data, 1.3, yKey-20, 'population', align=0.5, color=colSSO cgColorBar, nColors=256, range=resRang*1000, position=[0.92, 0.1, 0.98, 0.95], /Vertical, palette=palette cgText, /Normal, 0.86, 0.525, 'Residuals (mmag)', orient=90, align=0.5 cgPS_close, /png, Delete_PS=0 ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- VI -- Figures 2d Distro of Results -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; cgPS_open, xSize=22, ySize=26.5, Filename=dirPHOT+'euclid-p_vs_a-'+string(sigNoise,format='(F5.3)')+'-sig.eps', $ /Metric, /Decomposed, /Landscape, /Encapsulated, Language_Level=2, /Quiet ;--VI.1-- Parameters --------------------------------------------------------------------------; xBot = 0.15 xMid = 0.8 xTop = 0.9 yBot = 0.1 yTop = 0.95 pArr1 = [ xBot, yBot, xMid, yTop ] pArr2 = [ xMid, yBot, xTop, yTop ] ; cgLoadCT, 28, NCOLORS=256, RGB_TABLE=palette, /Brewer cgLoadCT, 39, NCOLORS=256, RGB_TABLE=palette ;--VI.2-- Uncertainty---------------------------------------------------------------------------; print, sigRang pp = where( uncP_vs_A[*,*,0,0] ge sigRang[0] and uncP_vs_A[*,*,0,0] le sigRang[1] ,nn ) print, 100.*nn/n_elements(uncP_vs_A[*,*,0,0]) loc =[0.001,0.5] uncImg = cgImgScl( uncP_vs_A[*,*,0,0], minValue=loc[0], maxValue=loc[1], stretch=6) cgImage, /NoErase, uncImg, /Normal, position=pArr1, palette=palette lvlArr = [0.00015, 0.00283, 0.0056] contour, pdsPDF, min_value=min(pdsPDF), max_value=max(pdsPDF), $ levels=lvlArr, Path_Info=path, Path_Xy=XY, /Closed, /Isotropic, /Path_Data_Coords nbCurve=n_elements(path) ; for kP=0,nbCurve-1 do begin yShi=[10, 3, 0.5] for kP=0,1 do begin pathIndex= [indgen(path(kP).N), 0] c = xy(*,path(kP).offset + pathIndex ) cgPlot, /OverPlot, c[0,*], c[1,*], color=colSSO ww= where( pdsPDF ge lvlArr[kP]) frac = string(total(pdsPDF[ww])*100.,format='(I3)') print, 'PDS : '+string(lvlArr[kP],format='(F5.3)')+' : '+frac+' %' yL = max(c[1,*],pMax) xL = c[0,pMax] cgText, /Data, xL, yL, frac+'%', color=colSSO endfor cgPlot, 0,vecP/3600., /NoData, /NoErase, /Normal, $ position=pArr1, $ xTickLen=0.02, yTickLen=0.04, $ xStyle=1, xRange=[amp.min, amp.max], $ yStyle=1, yRange=[period.min, period.max]/3600., /yLog, $ xTitle='Lightcurve amplitude (mag)', $ yTitle='Rotation period (h)' xKey=[0.8,1.1] yKey=130*[1,1] cgPlot, xKey, yKey, /OverPlot, color=colSSO cgText, /Data, 1.3, yKey+5, 'Known' , align=0.5, color=colSSO cgText, /Data, 1.3, yKey-20, 'population', align=0.5, color=colSSO cgColorBar, nColors=256, range=loc*1e3, position=[0.92, 0.1, 0.98, 0.95], /Vertical, palette=palette, /yLog cgText, /Normal, 0.86, 0.525, 'Dispersion (mmag)', orient=90, align=0.5 ; for k=0, period.N-1 do $ ; cgtext, 0.1+0.025*k, vecPout[k], strtrim(string(vecPout[k],format='(F5.1)'),2) cgPS_close, /png, Delete_PS=0 stop ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; ;--- TAG --- VII -- Illustration of LC in Euclid -----------------------; ;-----------------------------------------------------------------------------------------------; ;-----------------------------------------------------------------------------------------------; cgPS_open, xSize=16, ySize=26.5, Filename=dirPHOT+'euclid-simulc.eps', $ /Metric, /Decomposed, /Landscape, /Encapsulated, Language_Level=2, /Quiet kSSO = nbSSO-1 cgPlot, 0,0, /NoData, /Normal, $ xTitle='Time (min)', $ yTitle='Relative magnitude (arbitrarily shifted for clarity)', $ position=[0.15, 0.1, 0.99,0.99], $ xRange=[0,115], $ yRange=[0, 3.5] for kEx=0,5 do begin case kEx of 0: begin pMin = 2 pMax = 2.1 aMin = 0.13 aMax = 0.17 yMag = 0.3 end 1: begin pMin = 4 pMax = 4.1 aMin = 0.10 aMax = 0.15 yMag = .8 end 2: begin pMin = 4 pMax = 4.1 aMin = 0.70 aMax = 0.75 yMag = 1.5 end 3: begin pMin = 6 pMax = 6.1 aMin = 0.40 aMax = 0.45 yMag = 1.7 end 4: begin pMin = 10 pMax = 10.1 aMin = 0.20 aMax = 0.25 yMag = 2.4 end 5: begin pMin = 50 pMax = 55 aMin = 0.30 aMax = 0.35 yMag = 2.8 end else: endcase ; yMag=0 sel = where( stat.period/3600. gt pMin and stat.period/3600. lt pMax and $ stat.amp gt aMin and stat.amp lt aMax, nbSel ) xy=array_indices(stat, sel[0]) curLC = reform( bigLC[xy[0],*,xy[1]] ) curTime = ( obsStart+obsDIT/2. ) /60. curMag = reform( mag[xy[0],xy[1],*] ) curFilt = filtId print, kEx, nbSel, mean(curLC), stat[sel[0]].amp, stat[sel[0]].period/3600 cgPlot, /OverPlot, Time/60, yMag+curLC, color='Dark Gray' ; cgText, 70, yMag+curLC[-1], strTrim(string(kEx),2) cgText, 70, yMag+curLC[-1]+0.05, cgGreek('Delta')+'m: '+strTrim(string(stat[sel[0]].amp,format='(F3.1)'),2)+$ ', P:'+strTrim(string(stat[sel[0]].period/3600,format='(F4.1)'),2)+'h', $ charSize=1.5 for kP=0, 15 do begin case curFilt[kP] of 'VIS': begin locColor = 0. col='Black' sym='Filled Circle' end 'Y': begin ; locColor = 0.05 locColor = VmY[xy[0]] col='Green' sym='Open Square' end 'J': begin ; locColor = 0.1 locColor = VmJ[xy[0]] col='Cornflower Blue' sym='Open Circle' end 'H': begin ; locColor = 0.15 locColor = VmH[xy[0]] col='Red' sym='OpenUpTriangle' end endcase cgPlot, /OverPlot, Time/60, yMag+curLC-locColor/10., color=col cgPlot, /OverPlot, Time/60, yMag+curLC-locColor/10., color=col cgPlot, /OverPlot, Time/60, yMag+curLC-locColor/10., color=col range = where( time ge obsStart[kP] and time le (obsStart[kP]+obsDIT[kP])) loCmag = mean(curLC[range]) cgPlot, /OverPlot, curTime[kP], yMag+locMag-locColor/10, psym=sym, color=col ;, $ ; cgPlot, /OverPlot, curTime[kP], yMag+curMag[kP]+locColor, psym=sym, color=col;, $ ; Err_yLow=sigNoise, Err_yHigh=sigNoise endfor endfor xKey = [15, 15, 45, 45] ; yKey = [3.35, 3.25, 3.35, 3.25] yKey = [3.3, 3.2, 3.3, 3.2] yTxt = -0.035 xLen = 5 xShi = 7 for kF=0, 3 do begin cgPlot, /OverPlot, xKey[kF], yKey[kF], psym=symFilt[kF] , color=colFilt[kF] cgPlot, /OverPlot, xKey[kF]+[-1,1]*xLen, yKey[[kF,kF]], color=colFilt[kF] cgText, xKey[kF]+xShi, yKey[kF] + yTxt, listFilt[kF], color=colFilt[kF] endfor cgPS_close, /png, Delete_PS=0 end