;-Global settings root = '/observ/k2/' campaign = 'c1' doPlot = 1 ;-- generate plots of LC SetDecomposedState, 1, CurrentState=theState thisDevice = !D.Name ;-Selection criteria minFlux = 0. maxSide = 3. minSNR = 3. maxDEV = 3. ;-I/O Directories dirLC = root+campaign+'/LC/extraction/' dirEX = root+campaign+'/LC/single/' ;-Retrieve the target list & Loop over targets spawn, 'ls -1 '+dirLC+'*.txt', list nbAst = n_elements(list) ;forprint, indgen(nbAst), ' '+list, textout=2, subset=indgen(200)+1100 ;stop cKMperS = 299792.458D ;-km/s AUtoKM = 149598000.D ;-km cAUperS = cKMperS / AUtoKM ;-au/s cAUperD = cAUperS * 86400.D ;au/day ; for kAst=0, nbAst-1 do begin ; for kAst=1784, nbAst-1 do begin ;-c1-Sinon ; for kAst=1257, nbAst-1 do begin ;-c1-Antilochios for kAst=0, nbAst-1 do begin ;-c1-Amphilochios ;-------------------------------------------------------------------------------- ;-Identify target K2Name = file_basename( list(kAst) ) ssoName = strmid( K2Name, 0, strlen(K2Name)-7 ) dirName = ssoName underscore = strpos( ssoName, "_" ) if underscore ne -1 then strput, ssoName, ' ', underscore ssoInfo = voSsoDNet_resolver( ssoName ) ssoInfo = ssoInfo(0) if ssoInfo.number eq 0 or ssoInfo.Number ge 10000 then goto, nextSSO if ssoInfo.number eq 0 then begin dirAST = dirEX+dirName+'/' print, string(kAst,format='(I5)')+'/'+string(nbAst,format='(I5)')+': '+dirName ssoLabel = ' '+string(dirName,format='(A-20)')+' ' id = ssoInfo.name endif else begin dirAST = dirEX+string(ssoInfo.number,format='(I06)')+'-'+dirName+'/' print, string(kAst,format='(I5)')+'/'+string(nbAst,format='(I5)')+': ('+$ strtrim(string(ssoInfo.number,format='(I6)'),2)+') '+dirName ssoLabel = string(ssoInfo.number,format='(I6)')+' '+string(dirName,format='(A-20)') id = ssoInfo.number endelse if not file_exist(dirAst) then file_mkdir, dirAST ;-------------------------------------------------------------------------------- ;-Read Current Asteroid Lightcurve readcol, list(kAst), BJD, flux, dFlux, EPIC, flagSide, flagVar, flagSat, format='(D,F,F,L,F,I,I)', /silent nbPointIn = n_elements(EPIC) ;-------------------------------------------------------------------------------- ;-Extract GOOD photometry sel = where( flagVar ne 1 and $ ;-No variable star flagSat ne 1 and $ ;-No saturated flux ge minFlux and $ ;-Minimum flux flagSide le maxSide and $ ;-No SSO on the edges flux/dFlux ge minSNR, nSel ) ;-Minimum SNR ;-------------------------------------------------------------------------------- ;-Stat on rejection bad={all:0, sel:0, flux:0, snr:0, side:0, var:0, sat:0, dev:0} bad.all = nbPointIn bad.sel = nSel bad.flux = n_elements( where( flux le minFlux ) ) bad.side = n_elements( where( flagSide ge maxSide ) ) bad.snr = n_elements( where( flux/dFlux le minSNR ) ) bad.var = n_elements( where( flagVar eq 1 )) bad.sat = n_elements( where( flagSat eq 1 )) if nSel gt 1 then begin ;-------------------------------------------------------------------------------- ;-Convert BJD into UTC jd = dblarr(nbPointIn) forprint, BJD[sel], format='(D16.8)', textout='/tmp/dates.bjd',/NoComment, /Silent radec = voMiriade_callEphemCC( id, '/tmp/dates.bjd', obs='@Kepler', /web, /dateList, /JULIAN, tcoor=1, rPlane=1, mime='text',dump='/tmp/radec.eph' ) for kSel=0, nSel-1 do begin utc = bjd2utc( bjd[sel[kSel]], radec[kSel].ra.dec, radec[kSel].dec.dec) jd[sel[kSel]] = utc endfor ;------Get Geometry LC = generateLC(id, JD[sel], LTC=1, obs='@Kepler' ) ;------REDUCED MAGNITUDE -- DISTANCE dObs = sqrt(total( lc.obs^2, 2)) dSun = sqrt(total( lc.sun^2, 2)) corD = dObs*dOBs*dSun*dSun ;----- REDUCED MAG -- PHASE ANGLE dot = TOTAL( lc.OBs * lc.Sun, 2 ) phase = acos( dot/(dObs*dSun) ) 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 ) corP = (1-G)*phi1 + G*phi2 ;------Flux corrected from distance and Phase function flux11_ = flux[sel] * corD / mean(corD) flux110 = flux[sel] * corD * corP / mean(corD) / mean(corP) ;------Remove outlyers n110=n_elements(flux110) x110=indgen(n110) a=regress( x110, flux110, const=b) line=a[0]*x110+b residuals = flux110-line sig = stddev(residuals) dev = residuals/sig valid = where( abs(dev) le maxDEV, nbValid, comp=bad, ncomp=nbBad ) xAll =x110 jdAll = LC.JD yAll = flux110 if nbBad ne 0 then begin sel = sel[valid] nSel= nbValid flux11_ = flux11_[valid] flux110 = flux110[valid] LC.NUM = nbValid LC.JD = LC.JD[valid] LC.Sun = LC.Sun[valid,*] LC.Obs = LC.Obs[valid,*] n110=n_elements(flux110) x110=indgen(n110) endif unc = dFlux[sel] ;--- number of FOVs uEpic = uniq( EPIC[sel], sort(epic[sel])) nbEpic = n_elements(uEpic) ;------Write LC on disk ; print, dirAst+strtrim(string(id),2)+'.lcg' nameLC = strtrim(string(id),2)+'.lcg' openW, idOut, dirAst+nameLC, /get_lun printf, idOut, '1' iso = date_conv(mean( (lc.JD) ), 'FITS') printf, idOut, strtrim(string(lc.Num,format='(I)'),2)+' 0 '+$ strtrim(string(mean( unc ),format='(F12.4)'),2)+' '+$ strmid(iso,0,10)+' '+strtrim('Kmag',2)+' '+strtrim('K2',2) for kP=0, lc.Num-1 do begin printf, idOut, lc.JD(kP), $ flux11_[kP], $ lc.sun(kP,0), $ lc.sun(kP,1), $ lc.sun(kP,2), $ lc.obs(kP,0), $ lc.obs(kP,1), $ lc.obs(kP,2), $ format='(D14.6,1x,F12.4,6(1x,E16.9))' endfor close, idOut free_lun, idOut ;-Make figure if doPlot eq 1 then begin cgPS_open, Filename=+dirAst+strtrim(string(id),2)+'-time.eps', $ /metric, /decomposed, /landscape, $ xsize=30, ysize=20, language_level=2, /quiet cgPlot, jdAll, yAll minJD = min(lc.jd) cgPlot, jdAll-minJD, yAll, psym='filled circle', /NORMAL, $ xtitle='Time from '+strtrim(string(minJD,format='(D18.6)'),2)+' (days)', $ ytitle='Flux (electrons)', $ position=[0.1,0.08,0.9,0.9] cgplot, /over, LC.JD-minJD, flux110, psym='filled circle', color='Cornflower blue' cgPS_close, /png, Delete_PS=0 cgPS_open, Filename=+dirAst+strtrim(string(id),2)+'-index.eps', $ /metric, /decomposed, /landscape, $ xsize=30, ysize=20, language_level=2, /quiet cgPlot, jdAll, yAll cgPlot, xAll, yAll, psym='filled circle', /NORMAL, $ xtitle='Measurement Index', $ ytitle='Flux (electrons)', $ position=[0.1,0.08,0.9,0.9] cgplot, x110, a[0]*x110+b , /over, /addcmd cgplot, x110, a[0]*x110+b-sig, /over, /addcmd, color='gray' cgplot, x110, a[0]*x110+b+sig, /over, /addcmd, color='gray' cgplot, x110, a[0]*x110+b-sig*2, /over, /addcmd, color='gray' cgplot, x110, a[0]*x110+b+sig*2, /over, /addcmd, color='gray' cgplot, x110, a[0]*x110+b-sig*3, /over, /addcmd, color='gray' cgplot, x110, a[0]*x110+b+sig*3, /over, /addcmd, color='gray' cgplot, /over, x110, flux110, psym='filled circle', color='Cornflower blue' cgPS_close, /png, Delete_PS=0 endif ;----- Summar file openW, idSum, dirAST+"summary.dat", /get_lun printf, idSum, string(dirAST,format='(A-45)')+' '+$ string(ssoInfo.class,format='(A-10)')+' '+$ ssoLabel+' '+$ string(nbEpic,format='(I3)')+' '+$ string(nbPointIn,format='(I4)')+' '+$ string(nSel,format='(I3)') close, idSum free_lun, idSum ;----- Inversion files spawn, 'cp -r /observ/k2/template_ell/* '+dirAst spawn, 'sed -i -e "s/NUMBER/'+strtrim(string(id),2)+'/" '+dirAst+'list_ell_1' spawn, 'wget -q "http://astro.troja.mff.cuni.cz/projects/dafeed/php.php?script=view_data_file&lowell=1&asteroid_number='+strtrim(string(id),2)+'&file=Lowell_01.lcg" '+$ '--user dfuser --password df2dm4au -O '+dirAst+'lowell.lcg' spawn, 'echo 2 > '+dirAst+'lcs/'+nameLC spawn, "awk 'NR>1' "+dirAst+nameLC+' >> '+dirAst+'lcs/'+nameLC spawn, "awk 'NR>1' "+dirAst+'lowell.lcg >> '+dirAst+'lcs/'+nameLC endif else begin openW, idSum, dirAST+"summary.dat", /get_lun printf, idSum, string(dirAST,format='(A-45)')+' '+$ string(ssoInfo.class,format='(A-10)')+' '+$ ssoLabel+' '+$ string(0,format='(I3)')+' '+$ string(nbPointIn,format='(I4)')+' '+$ string(0,format='(I3)') close, idSum free_lun, idSum endelse nextSSO: endfor set_plot, thisDevice SetDecomposedState, theState end