#!/usr/bin/env python import unittest import pandas import numpy from math import * import sys, aa, ROOT #this actually works but test data from the `pip` packets will be used (so they can be not available upon installation) #from km3net_testdata import data_path #filename = data_path("astro/antares_coordinate_systems_benchmark.csv") #print("filename: ",filename) import os path = os.environ['AADIR']+"/externals/km3net-testdata/km3net_testdata/data/" deg = 180/pi #eqgqal SLA # DR,DD dp J2000.0 RA,Dec # DL,DB dp galactic longitude and latitude L2,B2 #all inputs in radians #Here: in deg def eqtogal(eq_ra,eq_dec): gal = ROOT.sla.eqgal( eq_ra/deg, eq_dec/deg ) return gal.DL*deg, gal.DB*deg def galtoeq(gal_lon,gal_lat): eq = ROOT.sla.galeq( eq_ra/deg, eq_dec/deg ) return eq.DR*deg, eq.DD*deg #sla_rdplan_(NewMJD, np, longitude, latitude, ra_moon, dec_moon, diam); #inputs: NewMJD, np, longitude, latitude #np =3 Moon, 0 - Sun #output ra_moon, dec_moon, diam #struct rdplan #double RA;double DEC;double DIAM; #rdplan( const double DATE,const int NP,const double ELONG,const double PHI ) def celestialtoeq(date,time,objname,det): yr,month,day = map( int, date.split('-') ) hr,mm,s = map( float, time.split(":") ) sec = int(s) ns = int(1e9*(s-sec)) ts = ROOT.TTimeStamp( yr, month, day, int(hr), int(mm), sec, ns, True ); mjd = ts.AsJulianDate() - 2400000.5 objtype = 0 if objname == "Moon": objtype = 3 elif objname == "Sun": objtype = 0 else: print("ERROR:",objname,"is unknown") celobj = ROOT.sla.rdplan(mjd, objtype, det.longitude, det.latitude) return celobj.RA*deg, celobj.DEC*deg def eqtoeq2000(date,time,ra,dec,direct = True): yr,month,day = map( int, date.split('-') ) hr,mm,s = map( float, time.split(":") ) sec = int(s) ns = int(1e9*(s-sec)) ts = ROOT.TTimeStamp( yr, month, day, int(hr), int(mm), sec, ns, True ); eq = ROOT.EquatorialCoords(ra/deg,dec/deg) eqnew = eq.correct_to_j2000(ts,not direct) ranew = ROOT.wrap( eqnew.ra()) return ranew*deg, eqnew.dec_deg() def localfromtrack(phi,theta): r = ROOT.Vec() r.set_angles(theta/deg, phi/deg) r = -r #to point to source phi = ROOT.wrap( r.phi()) return phi*deg, r.theta()*deg def equatorial(date,time,det,phi,theta,J200flag): yr,month,day = map( int, date.split('-') ) hr,mm,s = map( float, time.split(":") ) sec = int(s) ns = int(1e9*(s-sec)) ts = ROOT.TTimeStamp( yr, month, day, int(hr), int(mm), sec, ns, True ); r = ROOT.Vec() r.set_angles(theta/deg, phi/deg) eq = ROOT.EquatorialCoords( det, r, ts, J200flag ) ra = ROOT.wrap( eq.ra()) return ra*deg, eq.dec_deg() def localfromequatorial(date,time,det,ra,dec): yr,month,day = map( int, date.split('-') ) hr,mm,s = map( float, time.split(":") ) sec = int(s) ns = int(1e9*(s-sec)) ts = ROOT.TTimeStamp( yr, month, day, int(hr), int(mm), sec, ns, True ); ra_h,ra_m,ra_s = map(float, ra.split(':') ) dec_deg,dec_min,dec_sec = map(float, dec.split(':') ) dec_sign = 1 if (dec_deg < 0): dec_sign = -1 eq = ROOT.EquatorialCoords( 0, 0) eq.set_ra(int(ra_h),int(ra_m),ra_s) eq.set_dec(int(numpy.fabs(dec_deg)), int(dec_min), dec_sec, dec_sign) r = eq.track_direction(det,ts) phi = ROOT.wrap( r.phi()) return phi*deg, r.theta()*deg class TestAntaresBenchmark(unittest.TestCase): def setUp(self): # Normally, the detector coordinates would come from the .detx file # But in this case, we want to replicate ANTARES results and we # have to set the coordinats ourselves. The function set_lonlat # also computes the utm grid and convergence angle member of Det. #From KM3NeT coordinate document: https://docs.google.com/document/d/0B6l8SNtndcwaZzZHSFhfZWxsUEU #In the case of the ANTARES detector (centre of the 12-Line detector is at 42o 47.935'N and 6o 09.942'E, March 2010 #which corresponds to UTM zone 32 with Easting = 268221.6 meters and Northing = 4742381.9 meters #note that for convenience all ANTARES detector coordinates are given with 260000 m and 4740000 m subtracted for the Easting and Northing values, respectively. #UTM zone 32 with Easting = 268221.6 meters and Northing = 4742381.9 meters from http://antares-wiki.in2p3.fr/doku.php?id=astro_coordinatetransformation_howto&s[]=benchmark#benchmarks_for_coordinate_calculations #det_test = ROOT.Det() #det_test.utm_zone = 32 #det_test.utm_ref_easting = 260000 #det_test.utm_ref_northing = 4740000 #det_test.set_lonlat_from_utm() #print("test_detector_lon_lat_deg:",det_test.longitude*deg,det_test.latitude*deg) #det_test2 = ROOT.Det() #det_test2.utm_zone = 32 #det_test2.utm_ref_easting = 268221.6 #det_test2.utm_ref_northing = 4742381.9 #det_test2.set_lonlat_from_utm() #print("test_detector_lon_lat_deg:",det_test2.longitude*deg,det_test2.latitude*deg) self.det = ROOT.Det() #self.det.set_lonlat( 6.1657/deg, 42.798917/deg ) self.det.set_lonlat( (6+9.942/60.)/deg, (42+47.935/60.)/deg ) #from http://antares-wiki.in2p3.fr/doku.php?id=astro_coordinatetransformation_howto&s[]=benchmark#benchmarks_for_coordinate_calculations #print("ANTARES") #print("detector_lon_lat_deg:",self.det.longitude*deg,self.det.latitude*deg) #print("detector_UTM: ", self.det.utm_zone,self.det.utm_ref_easting, self.det.utm_ref_northing, self.det.utm_ref_z) #<--- not all is filled actually self.atol = 0.02 def test_antares_coordinate_systems_benchmark(self): table = pandas.read_csv(path+"astro/antares_coordinate_systems_benchmark.csv",comment='#') #used to reorder the ANTARES table #table_new = table[['date','time','phi','theta','azimuth','zenith','RA','DEC','RA-J2000','DEC-J2000','gal_lon','gal_lat']] #table_new.to_csv("antares_coordinate_systems_benchmark.csv",sep=",") table[['az_this','zen_this']] = table.apply(lambda x: localfromtrack(x.phi, x.theta), axis=1, result_type='expand') table[['ra2000_deg_this','dec2000_deg_this']] = table.apply(lambda x: equatorial(x.date, x.time, self.det, x.phi, x.theta, True), axis=1, result_type='expand') table[['ra_deg_this','dec_deg_this']] = table.apply(lambda x: equatorial(x.date, x.time, self.det, x.phi, x.theta, False), axis=1, result_type='expand') table[['gal_lon_deg_this','gal_lat_deg_this']] = table.apply(lambda x: eqtogal(x.ra2000_deg_this, x.dec2000_deg_this), axis=1, result_type='expand') numpy.testing.assert_allclose(table['zenith'],table['zen_this'],rtol=0,atol=self.atol,verbose=True) numpy.testing.assert_allclose(table['azimuth'],table['az_this'],rtol=0,atol=self.atol,verbose=True) numpy.testing.assert_allclose(table['DEC-J2000'],table['dec2000_deg_this'],rtol=0,atol=self.atol,verbose=True) numpy.testing.assert_allclose(table['RA-J2000'],table['ra2000_deg_this'],rtol=0,atol=self.atol,verbose=True) numpy.testing.assert_allclose(table['DEC'],table['dec_deg_this'],rtol=0,atol=self.atol,verbose=True) numpy.testing.assert_allclose(table['RA'],table['ra_deg_this'],rtol=0,atol=self.atol,verbose=True) numpy.testing.assert_allclose(table['gal_lat'],table['gal_lat_deg_this'],rtol=0,atol=self.atol,verbose=True) numpy.testing.assert_allclose(table['gal_lon'],table['gal_lon_deg_this'],rtol=0,atol=self.atol,verbose=True) def test_antares_objects_benchmark(self): table = pandas.read_csv(path+"astro/antares_astro_objects_benchmark.csv",comment='#') #used to reorder the ANTARES table #table_new = table[['date','time','object','RA-J2000','DEC-J2000','azimuth','zenith','phi','theta']] #table_new.to_csv("antares_astro_objects_benchmark.csv",sep=",") table[['phi_this','theta_this']] = table.apply(lambda x: localfromequatorial(x.date, x.time, self.det, x['RA-J2000'], x['DEC-J2000']), axis=1, result_type='expand') table[['az_this','zen_this']] = table.apply(lambda x: localfromtrack(x.phi_this, x.theta_this), axis=1, result_type='expand') numpy.testing.assert_allclose(table['phi'],table['phi_this'],rtol=0,atol=self.atol,verbose=True) numpy.testing.assert_allclose(table['theta'],table['theta_this'],rtol=0,atol=self.atol,verbose=True) numpy.testing.assert_allclose(table['zenith'],table['zen_this'],rtol=0,atol=self.atol,verbose=True) numpy.testing.assert_allclose(table['azimuth'],table['az_this'],rtol=0,atol=self.atol,verbose=True) def test_antares_SunMoon_benchmark(self): table = pandas.read_csv(path+"astro/antares_moon_sun_position_benchmark.csv",comment='#') #used to reorder the ANTARES table #table_new = table[['date','time','object','azimuth','zenith','RA','DEC','RA-J2000','DEC-J2000']] #table_new.to_csv("antares_moon_sun_position_benchmark.csv") table[['ra_deg_this','dec_deg_this']] = table.apply(lambda x: celestialtoeq(x.date, x.time, x['object'], self.det), axis=1, result_type='expand') table[['ra2000_deg_this','dec2000_deg_this']] = table.apply(lambda x: eqtoeq2000(x.date, x.time,x.ra_deg_this,x.dec_deg_this), axis=1, result_type='expand') numpy.testing.assert_allclose(table['DEC-J2000'],table['dec2000_deg_this'],rtol=0,atol=self.atol,verbose=True) numpy.testing.assert_allclose(table['RA-J2000'],table['ra2000_deg_this'],rtol=0,atol=self.atol,verbose=True) numpy.testing.assert_allclose(table['DEC'],table['dec_deg_this'],rtol=0,atol=self.atol,verbose=True) numpy.testing.assert_allclose(table['RA'],table['ra_deg_this'],rtol=0,atol=self.atol,verbose=True) class TestARCABenchmark(unittest.TestCase): def setUp(self): detfile = path+"detx/KM3NeT_-00000001_20171212.detx" self.det = ROOT.Det(detfile) self.det.set_lonlat_from_utm() #print("detector_file:",detfile) #print("detector_lon_lat_deg:",self.det.longitude*deg,self.det.latitude*deg) #print("detector_UTM: ", self.det.utm_zone,self.det.utm_ref_easting, self.det.utm_ref_northing, self.det.utm_ref_z) self.atol = 0.02 def test_ARCA_coordinate_systems_benchmark(self): table = pandas.read_csv(path+"astro/ARCA_coordinate_systems_benchmark.csv",comment='#') table[['az_this','zen_this']] = table.apply(lambda x: localfromtrack(x.phi, x.theta), axis=1, result_type='expand') table[['ra2000_deg_this','dec2000_deg_this']] = table.apply(lambda x: equatorial(x.date, x.time, self.det, x.phi, x.theta, True), axis=1, result_type='expand') table[['ra_deg_this','dec_deg_this']] = table.apply(lambda x: equatorial(x.date, x.time, self.det, x.phi, x.theta, False), axis=1, result_type='expand') table[['gal_lon_deg_this','gal_lat_deg_this']] = table.apply(lambda x: eqtogal(x.ra2000_deg_this, x.dec2000_deg_this), axis=1, result_type='expand') #print(table) #this was used to save new tables #table_new = table[['date','time','phi','theta','az_this','zen_this','ra_deg_this','dec_deg_this','ra2000_deg_this','dec2000_deg_this','gal_lon_deg_this','gal_lat_deg_this']] #table_new.to_csv("ARCA_coordinate_systems_benchmark.csv",sep=",") numpy.testing.assert_allclose(table['zenith'],table['zen_this'],rtol=0,atol=self.atol,verbose=True) numpy.testing.assert_allclose(table['azimuth'],table['az_this'],rtol=0,atol=self.atol,verbose=True) numpy.testing.assert_allclose(table['DEC-J2000'],table['dec2000_deg_this'],rtol=0,atol=self.atol,verbose=True) numpy.testing.assert_allclose(table['RA-J2000'],table['ra2000_deg_this'],rtol=0,atol=self.atol,verbose=True) numpy.testing.assert_allclose(table['DEC'],table['dec_deg_this'],rtol=0,atol=self.atol,verbose=True) numpy.testing.assert_allclose(table['RA'],table['ra_deg_this'],rtol=0,atol=self.atol,verbose=True) numpy.testing.assert_allclose(table['gal_lat'],table['gal_lat_deg_this'],rtol=0,atol=self.atol,verbose=True) numpy.testing.assert_allclose(table['gal_lon'],table['gal_lon_deg_this'],rtol=0,atol=self.atol,verbose=True) def test_ARCA_objects_benchmark(self): table = pandas.read_csv(path+"astro/ARCA_astro_objects_benchmark.csv",comment='#') table[['phi_this','theta_this']] = table.apply(lambda x: localfromequatorial(x.date, x.time, self.det, x['RA-J2000'], x['DEC-J2000']), axis=1, result_type='expand') table[['az_this','zen_this']] = table.apply(lambda x: localfromtrack(x.phi_this, x.theta_this), axis=1, result_type='expand') #this was used to save new tables #table_new = table[['date','time','object','RA-J2000','DEC-J2000',"az_this","zen_this",'phi_this','theta_this']] #table_new.to_csv("ARCA_astro_objects_benchmark.csv",sep=",") numpy.testing.assert_allclose(table['phi'],table['phi_this'],rtol=0,atol=self.atol,verbose=True) numpy.testing.assert_allclose(table['theta'],table['theta_this'],rtol=0,atol=self.atol,verbose=True) numpy.testing.assert_allclose(table['zenith'],table['zen_this'],rtol=0,atol=self.atol,verbose=True) numpy.testing.assert_allclose(table['azimuth'],table['az_this'],rtol=0,atol=self.atol,verbose=True) def test_ARCA_SunMoon_benchmark(self): table = pandas.read_csv(path+"astro/ARCA_moon_sun_position_benchmark.csv",comment='#') table[['ra_deg_this','dec_deg_this']] = table.apply(lambda x: celestialtoeq(x.date, x.time, x['object'], self.det), axis=1, result_type='expand') table[['ra2000_deg_this','dec2000_deg_this']] = table.apply(lambda x: eqtoeq2000(x.date, x.time,x.ra_deg_this,x.dec_deg_this), axis=1, result_type='expand') #this was used to save new tables #table_new = table[['date','time','object','azimuth','zenith','ra_deg_this','dec_deg_this','ra2000_deg_this','dec2000_deg_this']] #table_new.rename(columns = {'ra_deg_this':'RA','dec_deg_this':'DEC','ra2000_deg_this':'RA-J2000','dec2000_deg_this':'DEC-J2000'}, inplace = True) #table_new.to_csv("ARCA_moon_sun_position_benchmark.csv") numpy.testing.assert_allclose(table['DEC-J2000'],table['dec2000_deg_this'],rtol=0,atol=self.atol,verbose=True) numpy.testing.assert_allclose(table['RA-J2000'],table['ra2000_deg_this'],rtol=0,atol=self.atol,verbose=True) numpy.testing.assert_allclose(table['DEC'],table['dec_deg_this'],rtol=0,atol=self.atol,verbose=True) numpy.testing.assert_allclose(table['RA'],table['ra_deg_this'],rtol=0,atol=self.atol,verbose=True) class TestORCABenchmark(unittest.TestCase): def setUp(self): detfile = path+"detx/D_ORCA006_t.A02181836.p.A02181837.r.A02182001.detx" self.det = ROOT.Det(detfile) self.det.set_lonlat_from_utm() #print("detector_file:",detfile) #print("detector_lon_lat_deg:",self.det.longitude*deg,self.det.latitude*deg) #print("detector_UTM: ", self.det.utm_zone,self.det.utm_ref_easting, self.det.utm_ref_northing, self.det.utm_ref_z) self.atol = 0.02 def test_ORCA_coordinate_systems_benchmark(self): table = pandas.read_csv(path+"astro/ORCA_coordinate_systems_benchmark.csv",comment='#') table[['az_this','zen_this']] = table.apply(lambda x: localfromtrack(x.phi, x.theta), axis=1, result_type='expand') table[['ra2000_deg_this','dec2000_deg_this']] = table.apply(lambda x: equatorial(x.date, x.time, self.det, x.phi, x.theta, True), axis=1, result_type='expand') table[['ra_deg_this','dec_deg_this']] = table.apply(lambda x: equatorial(x.date, x.time, self.det, x.phi, x.theta, False), axis=1, result_type='expand') table[['gal_lon_deg_this','gal_lat_deg_this']] = table.apply(lambda x: eqtogal(x.ra2000_deg_this, x.dec2000_deg_this), axis=1, result_type='expand') #print(table) #this was used to save new tables #table_new = table[['date','time','phi','theta','az_this','zen_this','ra_deg_this','dec_deg_this','ra2000_deg_this','dec2000_deg_this','gal_lon_deg_this','gal_lat_deg_this']] #table_new.to_csv("ORCA_coordinate_systems_benchmark.csv",sep=",") numpy.testing.assert_allclose(table['zenith'],table['zen_this'],rtol=0,atol=self.atol,verbose=True) numpy.testing.assert_allclose(table['azimuth'],table['az_this'],rtol=0,atol=self.atol,verbose=True) numpy.testing.assert_allclose(table['DEC-J2000'],table['dec2000_deg_this'],rtol=0,atol=self.atol,verbose=True) numpy.testing.assert_allclose(table['RA-J2000'],table['ra2000_deg_this'],rtol=0,atol=self.atol,verbose=True) numpy.testing.assert_allclose(table['DEC'],table['dec_deg_this'],rtol=0,atol=self.atol,verbose=True) numpy.testing.assert_allclose(table['RA'],table['ra_deg_this'],rtol=0,atol=self.atol,verbose=True) numpy.testing.assert_allclose(table['gal_lat'],table['gal_lat_deg_this'],rtol=0,atol=self.atol,verbose=True) numpy.testing.assert_allclose(table['gal_lon'],table['gal_lon_deg_this'],rtol=0,atol=self.atol,verbose=True) def test_ORCA_objects_benchmark(self): table = pandas.read_csv(path+"astro/ORCA_astro_objects_benchmark.csv",comment='#') table[['phi_this','theta_this']] = table.apply(lambda x: localfromequatorial(x.date, x.time, self.det, x['RA-J2000'], x['DEC-J2000']), axis=1, result_type='expand') table[['az_this','zen_this']] = table.apply(lambda x: localfromtrack(x.phi_this, x.theta_this), axis=1, result_type='expand') #this was used to save new tables #table_new = table[['date','time','object','RA-J2000','DEC-J2000',"az_this","zen_this",'phi_this','theta_this']] #table_new.to_csv("ORCA_astro_objects_benchmark.csv",sep=",") numpy.testing.assert_allclose(table['phi'],table['phi_this'],rtol=0,atol=self.atol,verbose=True) numpy.testing.assert_allclose(table['theta'],table['theta_this'],rtol=0,atol=self.atol,verbose=True) numpy.testing.assert_allclose(table['zenith'],table['zen_this'],rtol=0,atol=self.atol,verbose=True) numpy.testing.assert_allclose(table['azimuth'],table['az_this'],rtol=0,atol=self.atol,verbose=True) def test_ORCA_SunMoon_benchmark(self): table = pandas.read_csv(path+"astro/ORCA_moon_sun_position_benchmark.csv",comment='#') table[['ra_deg_this','dec_deg_this']] = table.apply(lambda x: celestialtoeq(x.date, x.time, x['object'], self.det), axis=1, result_type='expand') table[['ra2000_deg_this','dec2000_deg_this']] = table.apply(lambda x: eqtoeq2000(x.date, x.time,x.ra_deg_this,x.dec_deg_this), axis=1, result_type='expand') #this was used to save new tables #table_new = table[['date','time','object','azimuth','zenith','ra_deg_this','dec_deg_this','ra2000_deg_this','dec2000_deg_this']] #table_new.rename(columns = {'ra_deg_this':'RA','dec_deg_this':'DEC','ra2000_deg_this':'RA-J2000','dec2000_deg_this':'DEC-J2000'}, inplace = True) #table_new.to_csv("ORCA_moon_sun_position_benchmark.csv") numpy.testing.assert_allclose(table['DEC-J2000'],table['dec2000_deg_this'],rtol=0,atol=self.atol,verbose=True) numpy.testing.assert_allclose(table['RA-J2000'],table['ra2000_deg_this'],rtol=0,atol=self.atol,verbose=True) numpy.testing.assert_allclose(table['DEC'],table['dec_deg_this'],rtol=0,atol=self.atol,verbose=True) numpy.testing.assert_allclose(table['RA'],table['ra_deg_this'],rtol=0,atol=self.atol,verbose=True) if __name__ == '__main__': unittest.main()