# You may want to reconfigure this according to your envrionment. 
STILTS?=stilts

# Define Healpix level

HEALPIX_LEVEL=6
PNG_X_RES=640
PNG_Y_RES=400

.PRECIOUS: %.systematic.vot %.average.vot

# Getting healpix for hsoy

TAPQUERY=$(STILTS) tapquery \
  tapurl='http://dc.g-vo.org/tap' \
  maxrec=200000000 \
  omode=out \
  out=$@ \
  ofmt=vot \
  executionduration=14400 

HEALPIXPLOT=$(STILTS) plot2sky \
		auxmap=inferno \
		auxlabel='$*'\
		auxvisible=true \
		legend=false \
		omode=out \
		ofmt=png \
		viewsys=galactic \
		out=$@ \
		xpix=$(PNG_X_RES) ypix=$(PNG_Y_RES) \
		projection=aitoff \
		layer_1=healpix \
			datalevel_1=$(HEALPIX_LEVEL)\
			datasys_1=equatorial \
			degrade_1=0 \
			combiner_1=sum \
			transparency_1=0 \
			healpix_1=hpx \
			in_1=$< \
			ifmt_1=votable \
			istream_1=true


%.systematic.vot:
	$(TAPQUERY) "SELECT \
    COUNT(*) as n, \
    ivo_healpix_index ($(HEALPIX_LEVEL), hsoy.raj2000, hsoy.dej2000) AS hpx, \
    AVG(hsoy.$*-ppmxl.$*) AS obs \
    FROM hsoy.main AS hsoy \
    JOIN ppmxl.main AS ppmxl \
    USING (ipix)\
    GROUP BY hpx\
    "

%.average.vot:
	$(TAPQUERY) adql="SELECT \
    COUNT(*) as n, \
    ivo_healpix_index ($(HEALPIX_LEVEL), raj2000, dej2000) AS hpx, \
    AVG($*) AS obs \
    FROM hsoy.main AS hsoy \
    GROUP BY hpx\
    "

%.png: %.vot
	$(HEALPIXPLOT) value_1='obs'


all: autocorr_distances.vot e_pmde.average.png e_pmra.average.png \
	gmag_dist.png e_pm_vs_gmag.png e_pos_vs_gmag.png densityplot.png \
	pair_density.png

densityplot.png: e_pmde.average.vot
	$(HEALPIXPLOT) value_1='n'

e_pm_vs_gmag.vot:
	$(TAPQUERY) adql="select in_unit(avg(e_pmra), 'mas/yr') as e_pmra_avg, \
		in_unit(avg(e_pmde), 'mas/yr') as e_pmde_avg, \
		round(phot_g_mean_mag) as bin \
		from hsoy.main \
		group by bin"

e_pos_vs_gmag.vot:
	$(TAPQUERY) adql="select in_unit(avg(e_deepDE), 'mas') as e_ra_avg, \
		in_unit(avg(e_raepRA), 'mas') as e_dec_avg, \
		round(phot_g_mean_mag) as bin \
		from hsoy.main \
		group by bin"


autocorr_distances.vot:
	$(TAPQUERY) adql="select \
		in_unit(distance(point('', a.raj2000, a.dej2000), \
			point('', b.raj2000, b.dej2000)), 'arcsec') as dist \
		from hsoy.main as a \
		join (select all top 1000000000 * from hsoy.main) as b \
		on 1=contains(point('', a.raj2000, a.dej2000), \
				circle('', b.raj2000, b.dej2000, 5./3600))\
			and (a.ipix!=b.ipix and a.comp!=b.comp)"

autocorr_gaia_distances.vot:
	$(TAPQUERY) adql="select \
		in_unit(distance(point('', a.raj2000, a.dej2000), \
			point('', b.raj2000, b.dej2000)), 'arcsec') as dist \
		from gaia.dr1 as a \
		join (select all top 100000000 * from gaia.dr1) as b \
		on 1=contains(point('', a.raj2000, a.dej2000), \
				circle('', b.raj2000, b.dej2000, 5./3600))\
			and a.source_id!=b.source_id"

gmag_dist.vot:
	$(TAPQUERY) adql="select \
		round(phot_g_mean_mag*10)/10 as bin, \
		count(*) as n \
	from hsoy.main \
	group by bin"

b2mag_dist_ppmxl.vot:
	$(TAPQUERY) adql="select \
		round(b2mag) as bin, \
		count(*) as n \
	from ppmxl.main \
	group by bin"

gmag_dist_gaia.vot:
	$(TAPQUERY) adql="select \
		round(phot_g_mean_mag*10)/10 as bin, \
		count(*) as n \
	from gaia.dr1 \
	group by bin"

pair_density.vot:
	$(TAPQUERY) adql="SELECT \
    COUNT(*) as n, \
    ivo_healpix_index ($(HEALPIX_LEVEL), raj2000, dej2000) AS hpx, \
    (SUM(comp)+SUM(clone))/(1.0*COUNT(*)) as obs \
    FROM hsoy.main AS hsoy \
    GROUP BY hpx\
    "

tgas_comparison.vot:
	$(TAPQUERY) adql="SELECT \
    COUNT(*) as n, \
    ivo_healpix_index ($(HEALPIX_LEVEL), hsoy.raj2000, hsoy.dej2000) AS hpx, \
    SUM(tgas.pmra-hsoy.pmra*3.6e6)/count(*) as diff_pmra,\
    SUM(tgas.pmdec-hsoy.pmde*3.6e6)/count(*) as diff_pmde,\
    SQRT(SUM(POWER(tgas.pmra-hsoy.pmra*3.6e6, 2)))/(count(*)+1.) as scatter_pmra,\
    SQRT(SUM(POWER(tgas.pmdec-hsoy.pmde*3.6e6, 2)))/(count(*)+1.) as scatter_pmde\
    FROM tgas.main AS tgas \
    JOIN hsoy.main AS hsoy \
    ON (hsoy.gaia_id=tgas.source_id) \
    GROUP BY hpx\
    "

ucac_radiffs.vot:
	$(TAPQUERY) adql="SELECT \
    in_unit(avg(uc.pmra-hsoy.pmra), 'mas/yr') AS pmradiff, \
    in_unit(avg(uc.pmde-hsoy.pmde), 'mas/yr') AS pmdediff, \
    count(*) as n, \
    ivo_healpix_index (6, raj2000, dej2000) AS hpx \
    FROM hsoy.main AS hsoy \
    JOIN ucac5.main as uc \
    ON (uc.source_id=hsoy.gaia_id) \
    WHERE comp IS NULL \
    AND clone IS NULL \
    GROUP BY hpx\
    "

comp_density.vot:
	$(TAPQUERY) adql="SELECT \
    COUNT(*) AS obs, \
    ivo_healpix_index ($(HEALPIX_LEVEL), raj2000, dej2000) AS hpx \
    FROM hsoy.main AS hsoy \
    WHERE comp IS NOT NULL \
    GROUP BY hpx\
    "

highpm_density.vot:
	$(TAPQUERY) adql="SELECT \
    COUNT(*) AS obs, \
    ivo_healpix_index ($(HEALPIX_LEVEL), raj2000, dej2000) AS hpx \
    FROM hsoy.main AS hsoy \
    WHERE sqrt(pmra*pmra+pmde*pmde)>150/3.6e6 \
    	AND no_sc IS NULL  and comp IS NULL and clone IS NULL \
    GROUP BY hpx\
    "

pm_histogram.vot:
	$(TAPQUERY) adql="SELECT \
		round(pmtot/10)*10 as bin, \
		count(*) as n \
		FROM ( \
			SELECT sqrt(pmra*pmra+pmde*pmde)*3.6e6 as pmtot \
			FROM hsoy.main) AS q \
		group by bin"

pm_histogram_cleaned.vot:
	$(TAPQUERY) adql="SELECT \
		round(pmtot/10)*10 as bin, \
		count(*) as n \
		FROM ( \
			SELECT sqrt(pmra*pmra+pmde*pmde)*3.6e6 as pmtot \
			FROM hsoy.main\
			WHERE no_sc IS NULL and comp IS NULL and clone IS NULL) AS q \
		group by bin"

pm_histograms.png: pm_histogram.vot pm_histogram_cleaned.vot
	$(STILTS) plot2plane legend=false omode=out ofmt=png out=$@ \
		title="All-sky" \
		xpix=$(PNG_Y_RES) ypix=$(PNG_Y_RES) \
		ylog=True xlog=True\
		xlabel="PM bin [mas/yr, bin size=10]" \
		xmax=4000 \
	layer1=mark \
		color1=blue \
		in1=pm_histogram.vot \
		x1=bin \
		y1=n \
	layer2=mark \
		in2=pm_histogram_cleaned.vot \
		x2=bin \
		y2=n

pm_histogram_north.vot:
	$(TAPQUERY) adql="SELECT \
		round(pmtot/10)*10 as bin, \
		count(*) as n \
		FROM ( \
			SELECT sqrt(pmra*pmra+pmde*pmde)*3.6e6 as pmtot \
			FROM hsoy.main \
			WHERE dej2000>0) AS q \
		group by bin"

pm_histogram_cleaned_north.vot:
	$(TAPQUERY) adql="SELECT \
		round(pmtot/10)*10 as bin, \
		count(*) as n \
		FROM ( \
			SELECT sqrt(pmra*pmra+pmde*pmde)*3.6e6 as pmtot \
			FROM hsoy.main\
			WHERE no_sc IS NULL and comp IS NULL and clone IS NULL AND dej2000>0) AS q \
		group by bin"

pm_histograms_north.png: pm_histogram_north.vot pm_histogram_cleaned_north.vot
	$(STILTS) plot2plane legend=false omode=out ofmt=png out=$@ \
		title="Northern sky" \
		xpix=$(PNG_Y_RES) ypix=$(PNG_Y_RES) \
		ylog=True xlog=True\
		xlabel="PM bin [mas/yr, bin size=10]" \
		xmax=4000 \
	layer1=mark \
		color1=blue \
		in1=pm_histogram_north.vot \
		x1=bin \
		y1=n \
	layer2=mark \
		in2=pm_histogram_cleaned_north.vot \
		x2=bin \
		y2=n

incredibly_fast.vot: Makefile
	$(TAPQUERY) adql="SELECT h.ipix, h.comp, \
			sqrt(p.pmra*p.pmra+p.pmde*p.pmde) as ppmtot,  \
			sqrt(s.pmra*s.pmra+s.pmde*s.pmde) as spmtot,  \
			atan2(p.pmde, p.pmra) as pangle,  \
			atan2(s.pmra, s.pmde) as sangle \
FROM \
(SELECT TOP 50000 \
		ipix, comp, pmra, pmde \
		FROM hsoy.main \
		WHERE pmra>0.75/3600 or pmde>0.75/3600) as h \
JOIN ppmxl.main as p \
USING (ipix) \
LEFT OUTER JOIN supercosmos.sources as s \
		ON 1=CONTAINS(POINT('ICRS', s.raj2000, s.dej2000), \
			CIRCLE('ICRS', p.raj2000, p.dej2000, 3/3600.))"
		
	
clean:
	rm -f *.vot *.png