fix input for c to python

StrayLight.cpp

@@ -144,7 +144,7 @@ public:
 	}
 };
 
-void Init();
+void Init(const char *deFn, const char *PSFFn, const char *RFn, const char *zolFn);
 double PST(double Azimuth,double Angle);//73*45，行是方位角角5度一分，列是夹角2度一分，PST每个数占16字符，每行720+2字符
 
 XYZ Spherical2Cartesian(double latitude,double longitude);//天球转三维，参数为弧度
@@ -155,10 +155,10 @@ XYZ Equatorial2Ecliptic(XYZ Equ);//赤道转黄道
 double ShelterEarth(XYZ sat,XYZ star);//卫星是否会被地球遮挡
 double ShelterPlate(XYZ sp,XYZ ob,XYZ py,double fai,double op,double sum);//亮星是否会被盖板遮挡
 
-void PointSource(double jd, double sat[3], double ob[3], double py[3], double E[7]);
+void PointSource(double jd, double sat[3], double ob[3], double py[3], double E[7], char* starMagFn);
 void EarthShine(double ju, double sat[3], double ob[3], double py[3], double E[7]);
 void Zodiacal(double ju, double ob[3], double E[7]);//琢磨半晌形参还是用赤道坐标系吧，内部再转换黄道！
-void Calculate(double ju, double sat[3], double ob[3], double py[3], double E[7]);
+void Calculate(double ju, double sat[3], double ob[3], double py[3], double E[7], char* starMagFn);
 
 void ComposeY(double ob[3], double py1[3], double py2[3]);
 void Watt2Photon(double E[7], double P[7]);//照度w/m2转光子数ph/(s*100um2)
@@ -169,11 +169,15 @@ const int Rrow = 180, Rcol = 288;
 double R[Rrow][Rcol];
 double ZL[37][19];
 
-string de405("DE405");
-Ephcom ephcom(de405);
-void Init()
+
+Ephcom ephcom = Ephcom();
+void Init(const char *deFn, const char *PSFFn, const char *RFn, const char *zolFn)
 {
-    ifstream infile("PST");
+//    string de405("DE405");
+    string de405(deFn);
+    ephcom.init(de405);
+//    ifstream infile("PST");
+    ifstream infile(PSFFn);
     string s;
     for (int i = 0; i < 73; i++)
     {
@@ -185,7 +189,8 @@ void Init()
         }
     }
     infile.close();
-    infile.open("R");
+//    infile.open("R");
+    infile.open(RFn);
     for (int i = 0; i < Rrow; i++)
     {
         for (int j = 0; j < Rcol; j++)
@@ -194,7 +199,8 @@ void Init()
         }
     }
     infile.close();
-    infile.open("Zodiacal");
+//    infile.open("Zodiacal");
+    infile.open(zolFn);
     for (int i = 0; i < 37; i++)
     {
         for (int j = 0; j < 19; j++)
@@ -288,7 +294,7 @@ double ShelterPlate(XYZ sp,XYZ ob,XYZ py,double fai,double op,double sum)//亮
     }
 }
 
-void PointSource(double ju, double sat[3], double ob[3], double py[3], double E[7])
+void PointSource(double ju, double sat[3], double ob[3], double py[3], double E[7], char* starMagFn)
 {
     XYZ SAT = XYZ(sat[0], sat[1], sat[2]);
     XYZ OB = XYZ(ob[0], ob[1], ob[2]);
@@ -392,7 +398,7 @@ void PointSource(double ju, double sat[3], double ob[3], double py[3], double E[
 
     string s;
     double Estar[7] = { 0, 0, 0, 0, 0, 0, 0 };
-    ifstream infile("BrightGaia_with_csst_mag");
+    ifstream infile(starMagFn);
     while (getline(infile, s))
     {
         stringstream ss(s);
@@ -674,10 +680,10 @@ void Zodiacal(double ju, double ob[3], double E[7])
     // cout << endl;
 }
 
-void Calculate(double ju, double sat[3], double ob[3], double py[3], double E[7])
+void Calculate(double ju, double sat[3], double ob[3], double py[3], double E[7], char* starMagFn)
 {
     double Ep[7],Ee[7],Ez[7];
-    PointSource(ju, sat,ob,py,Ep);
+    PointSource(ju, sat,ob,py,Ep, starMagFn);
     EarthShine(ju, sat,ob,py,Ee);
     Zodiacal(ju,ob,Ez);
     // cout << "总杂散光：";

ephcom_wrapper.cpp

 #include "ephcom_wrapper.hpp"
 
-Ephcom::Ephcom(std::string& DE_filename){
+Ephcom::Ephcom(){
     this->infp = NULL;
     this->ORIG = EPHCOM_EARTH;
-    this->init(DE_filename);
+//    this->init(DE_filename);
 }
 
 Ephcom::~Ephcom(){

ephcom_wrapper.hpp

@@ -43,7 +43,7 @@ private:
 	int ORIG;			// origin used to calculation positions of Sun, Moon and planets. Default: EPHCOM_EARTH
 public:
 
-	Ephcom(std::string& DE_filename);
+	Ephcom();
     ~Ephcom();
 
     void init(std::string& DE_filename);

straylight.py

@@ -15,6 +15,8 @@ bandRange = {'nuv':[2504.0,3230.0],'u':[3190.0,4039.0],'g':[3989.0,5498.0],'r':[
 Instrument_dir = '/Users/zhangxin/Work/SlitlessSim/CSST_SIM/CSST_C6/straylight/straylight/Instrument/'
 SpecOrder = ['-2','-1','0','1','2']
 
+filterMirrorEff = {'nuv':0.54,'u':0.68,'g':0.8,'r':0.8, 'i':0.8, 'z':0.8,'y':0.8}
+
 def transRaDec2D(ra, dec):
     x1 = np.cos(dec / 57.2957795) * np.cos(ra / 57.2957795);
     y1 = np.cos(dec / 57.2957795) * np.sin(ra / 57.2957795);
@@ -91,15 +93,22 @@ class StrayLight(object):
         self.ecliptic = self.equator.transform_to('barycentrictrueecliptic')
         self.pointing = transRaDec2D(radec[0], radec[1])
         self.slcdll=ctypes.CDLL('./libstraylight.dylib')
-        self.slcdll.Calculate.argtypes =[ctypes.c_double ,ctypes.POINTER(ctypes.c_double), ctypes.POINTER(ctypes.c_double),ctypes.POINTER(ctypes.c_double), ctypes.POINTER(ctypes.c_double)]
+        self.slcdll.Calculate.argtypes =[ctypes.c_double ,ctypes.POINTER(ctypes.c_double), ctypes.POINTER(ctypes.c_double),ctypes.POINTER(ctypes.c_double), ctypes.POINTER(ctypes.c_double), ctypes.c_char_p]
 
-        self.slcdll.PointSource.argtypes =[ctypes.c_double ,ctypes.POINTER(ctypes.c_double), ctypes.POINTER(ctypes.c_double),ctypes.POINTER(ctypes.c_double), ctypes.POINTER(ctypes.c_double)]
+        self.slcdll.PointSource.argtypes =[ctypes.c_double ,ctypes.POINTER(ctypes.c_double), ctypes.POINTER(ctypes.c_double),ctypes.POINTER(ctypes.c_double), ctypes.POINTER(ctypes.c_double), ctypes.c_char_p]
 
         self.slcdll.EarthShine.argtypes =[ctypes.c_double ,ctypes.POINTER(ctypes.c_double), ctypes.POINTER(ctypes.c_double),ctypes.POINTER(ctypes.c_double), ctypes.POINTER(ctypes.c_double)]
 
         self.slcdll.Zodiacal.argtypes =[ctypes.c_double ,ctypes.POINTER(ctypes.c_double), ctypes.POINTER(ctypes.c_double)]
         self.slcdll.ComposeY.argtypes=[ctypes.POINTER(ctypes.c_double),ctypes.POINTER(ctypes.c_double),ctypes.POINTER(ctypes.c_double)]
-        self.slcdll.Init()
+        self.slcdll.Init.argtypes=[ctypes.c_char_p,ctypes.c_char_p,ctypes.c_char_p,ctypes.c_char_p]
+        self.deFn = "DE405"
+        self.PSTFn = "PST"
+        self.RFn = "R"
+        self.ZolFn = "Zodiacal"
+        self.brightStarTabFn = "BrightGaia_with_csst_mag"
+
+        self.slcdll.Init(str.encode(self.deFn),str.encode(self.PSTFn),str.encode(self.RFn),str.encode(self.ZolFn))
 
     def getFilterAndCCD_Q(self, filter = 'i'):
         ccd_fn = Instrument_dir + 'ccd/' + filterCCD[filter] + '.txt'
@@ -232,7 +241,7 @@ class StrayLight(object):
 
 
         star_e1 = (ctypes.c_double*7)()
-        self.slcdll.PointSource(self.jtime,sat,ob,py,star_e1)
+        self.slcdll.PointSource(self.jtime,sat,ob,py,star_e1, str.encode(self.brightStarTabFn))
 
         band_star_e1 = star_e1[:][filterIndex[filter]]
 
@@ -240,56 +249,7 @@ class StrayLight(object):
         
         return pix_star_e1
 
-
-
-    def caculateStrayLightFilter(self, filter = 'i', pixel_size_phy = 10 ):
-        sat = (ctypes.c_double*3)()
-        sat[:] = self.sat
-        ob = (ctypes.c_double*3)()
-        ob[:]=self.pointing
-        
-        
-        py1 = (ctypes.c_double*3)()
-        py2 = (ctypes.c_double*3)()
-        self.slcdll.ComposeY(ob,py1,py2)
-
-
-        earth_e1 = (ctypes.c_double*7)()
-        self.slcdll.EarthShine(self.jtime,sat,ob,py1,earth_e1)
-        earth_e2 = (ctypes.c_double*7)()
-        self.slcdll.EarthShine(self.jtime,sat,ob,py2,earth_e2)
-        zodical_e = (ctypes.c_double*7)()
-        self.slcdll.Zodiacal(self.jtime,ob,zodical_e)
-        
-        band_earth_e1 = earth_e1[:][filterIndex[filter]]
-        band_earth_e2 = earth_e2[:][filterIndex[filter]]
-        # band_earth_e1 = 0
-        # band_earth_e2 = 0
-        band_zodical_e = zodical_e[:][filterIndex[filter]]
-        
-        q=self.getFilterAndCCD_Q(filter=filter)
-        p_lambda = filterPivotWave[filter]
-        c = cons.c.value
-        h = cons.h.value
-        pix_earth_e1 = band_earth_e1/(h*c/(p_lambda*1e-10))*pixel_size_phy*1e-6*pixel_size_phy*1e-6*q
-        pix_earth_e2 = band_earth_e2/(h*c/(p_lambda*1e-10))*pixel_size_phy*1e-6*pixel_size_phy*1e-6*q
-        pix_zodical_e = band_zodical_e/(h*c/(p_lambda*1e-10))*pixel_size_phy*1e-6*pixel_size_phy*1e-6*q
-
-        # star_e1 = (ctypes.c_double*7)()
-        # self.slcdll.PointSource(self.jtime,sat,ob,py1,star_e1)
-        # # star_e2 = (ctypes.c_double*7)()
-        # # self.slcdll.PointSource(self.jtime,sat,ob,py2,star_e2)
-        # band_star_e1 = star_e1[:][filterIndex[filter]]
-        # # band_star_e2 = star_e2[:][filterIndex[filter]]
-        # pix_star_e1 = band_star_e1/(h*c/(p_lambda*1e-10))*pixel_size_phy*1e-6*pixel_size_phy*1e-6*q
-
-
-        pix_star_e1 = 0
-        # pix_star_e2 = band_star_e2/(h*c/(p_lambda*1e-10))*pixel_size_phy*1e-6*pixel_size_phy*1e-6*q
-        
-        return pix_earth_e1+pix_zodical_e+pix_star_e1, pix_zodical_e, pix_earth_e1, pix_earth_e2
-
-    def caculateStrayLightGrating(self, grating = 'GU', pixel_size_phy = 10, normFilter = 'g',  aper = 2, pixelsize = 0.074):
+    def caculateEarthshineGrating(self, grating = 'GU', pixel_size_phy = 10, normFilter = 'g',  aper = 2, pixelsize = 0.074):
         sat = (ctypes.c_double*3)()
         sat[:] = self.sat
         ob = (ctypes.c_double*3)()
@@ -310,7 +270,12 @@ class StrayLight(object):
         
         band_earth_e1 = earth_e1[:][filterIndex[normFilter]]
         band_earth_e2 = earth_e2[:][filterIndex[normFilter]]
-        band_earth_e = np.min([band_earth_e1, band_earth_e2])
+        band_earth_e = band_earth_e2
+        py = py2[:]
+        if band_earth_e1<band_earth_e2:
+            band_earth_e = band_earth_e1
+            py = py1[:]
+        # band_earth_e = np.min([band_earth_e1, band_earth_e2])
 
         # band_earth_e1 = 0
         # band_earth_e2 = 0
@@ -331,7 +296,7 @@ class StrayLight(object):
 
         lamb_earth = earthshine_spec['WAVELENGTH']
         flux_earth = earthshine_spec['FLUX']*pix_earth_e/earthshine_v
-        print(pix_earth_e,earthshine_v)
+        # print(pix_earth_e,earthshine_v)
         earth_v_grating = 0
         for s_order in SpecOrder:
             thpFn = Instrument_dir + 'sls_conf/' + grating + '.Throughput.' + s_order + 'st.fits'
@@ -340,11 +305,64 @@ class StrayLight(object):
             thp = thpFn_i(lamb_earth)
             beamsEarth = np.trapz(flux_earth*thp,lamb_earth)* math.pi*aper*aper/4 * pixelsize * pixelsize
             earth_v_grating = earth_v_grating + beamsEarth
-            print(beamsEarth)
+            # print(beamsEarth)
 
         # print(earthshine_v, pix_earth_e, earth_v_grating)
-        return earth_v_grating
+        return earth_v_grating, py
+
+    def caculateStarLightGrating(self, grating = 'GU', pointYaxis = np.array([1,1,1]), pixel_size_phy = 10 ):
+        sat = (ctypes.c_double*3)()
+        sat[:] = self.sat
+        ob = (ctypes.c_double*3)()
+        ob[:]=self.pointing
+        
+        
+        py = (ctypes.c_double*3)()
+        py[:] = pointYaxis
+
+        # q=self.getFilterAndCCD_Q(filter=filter)
+        # p_lambda = filterPivotWave[filter]
+        c = cons.c.value
+        h = cons.h.value
+
+
+        star_e1 = (ctypes.c_double*7)()
+        self.slcdll.PointSource(self.jtime,sat,ob,py,star_e1, str.encode(self.brightStarTabFn))
+
+        filterPivotWaveList = np.zeros(7)
+        bandRangeList = np.zeros(7)
+        filterMirrorEffList = np.zeros(7)
+        filterNameList = list(filterPivotWave.keys())
+        for i in np.arange(7):
+            filterPivotWaveList[i] = filterPivotWave[filterNameList[i]]
+            filterMirrorEffList[i] = filterMirrorEff[filterNameList[i]]
+            brange = bandRange[filterNameList[i]]
+            bandRangeList[i] = brange[1] - brange[0]
+        filterFlux_lamb = star_e1[:]/bandRangeList/filterMirrorEffList/(h*c/(filterPivotWaveList*1e-10))
+        filterFlux_lambi = interpolate.interp1d(filterPivotWaveList,filterFlux_lamb,fill_value="extrapolate")
+
+        lamb_g = np.arange(bandRange[grating][0], bandRange[grating][1],1)
+        flux_g = filterFlux_lambi(lamb_g)
+        # flux_total_g = np.trapz(flux_g,lamb_g)
+        zoldialcal_grating = 0
+        for s_order in SpecOrder:
+            thpFn = Instrument_dir + 'sls_conf/' + grating + '.Throughput.' + s_order + 'st.fits'
+            thp_ = Table.read(thpFn)
+            thpFn_i = interpolate.interp1d(thp_['WAVELENGTH'], thp_['SENSITIVITY'])
+            thp = thpFn_i(lamb_g)
+            beamsZol = np.trapz(flux_g*thp,lamb_g)*pixel_size_phy*1e-6*pixel_size_phy*1e-6
+            zoldialcal_grating = zoldialcal_grating + beamsZol
+            # print(beamsZol)
+
+
+
+        # band_star_e1 = star_e1[:][filterIndex[filter]]
+
+        # pix_star_e1 = band_star_e1/(h*c/(p_lambda*1e-10))*pixel_size_phy*1e-6*pixel_size_phy*1e-6*q
         
+        return zoldialcal_grating
+
+
     
     def calculatSkylightBySpec(self, specType = 'earthshine', filter = 'g', aper = 2, pixelsize = 0.074, s = 2000, e = 11000):
         specFn = ''
@@ -396,7 +414,7 @@ class StrayLight(object):
 
         longitude = longitude - (sun_equtor.lon*u.degree).value
         longitude = np.abs(longitude)
-        print((sun_equtor.lon*u.degree).value)
+        # print((sun_equtor.lon*u.degree).value)
 
         if (longitude > 180):
             longitude = 360 - longitude
@@ -421,6 +439,20 @@ class StrayLight(object):
         # print("factor:", v_mag, lo, la)
 
         return Table(np.array([spec[:,0], spec[:,1]*f_ration]).T,names=('WAVELENGTH', 'FLUX')), v_mag
+    
+    def caculateStrayLightFilter(self, filter = 'i', pixel_size_phy = 10, pixel_scale = 0.074, sun_pos = np.array([0,0,0])):
+        e1,py = self.caculateEarthShineFilter(filter = filter, pixel_size_phy = pixel_size_phy)
+        e2, _ = self.calculateZaodiacalFilter2(filter = 'i', sun_pos=sun_pos, pixelsize = pixel_scale)
+        e3 = self.caculateStarLightFilter(filter = 'i',pointYaxis = py, pixel_size_phy = pixel_size_phy)
+        
+        return e1+e2+e3
+    
+    def caculateStrayLightGrating(self, grating = 'GI', pixel_size_phy = 10, normFilter_es = 'g', sun_pos = np.array([0,0,0])):
+        e1,py = self.caculateEarthshineGrating(grating = grating, pixel_size_phy = pixel_size_phy, normFilter = normFilter_es)
+        e2 = self.caculateStarLightGrating(grating = grating, pointYaxis = py)
+        spec, _ = self.calculateZaodicalSpec(longitude = self.ecliptic.lon.value, latitude = self.ecliptic.lat.value, sun_pos = sun_pos)
+        
+        return e1+e2, spec
  
 
 def testZodiacal(lon = 285.04312526255366, lat = 30.):
@@ -478,19 +510,26 @@ for i in np.arange(surveylist.shape[0]):
     sky_pix[i,4] = earthZenithAngle
     sky_pix[i,5] = isInSunSide
 
-    e1,py = sl.caculateEarthShineFilter(filter = 'i')
+    e1_,py = sl.caculateEarthShineFilter(filter = 'i')
     # e2, e2_ = sl.calculateZodiacalFilter1(filter = 'i')
     e3, v_mag = sl.calculateZaodiacalFilter2(filter = 'i', sun_pos=surveylist[i,6:9])
-    e4 = sl.caculateStarLightFilter(filter = 'i',pointYaxis = py)
+    # e4 = sl.caculateStarLightFilter(filter = 'i',pointYaxis = py)
     # e4 = 0 
 
-    # e4 = sl.caculateStrayLightGrating(grating = 'GI', pixel_size_phy = 10, normFilter = 'g')
+    e1,py = sl.caculateEarthshineGrating(grating = 'GI', pixel_size_phy = 10, normFilter = 'g')
+
+    # e2 = sl.caculateStarLightGrating(grating = 'GV', pointYaxis = py)
+    e2 = sl.caculateStarLightGrating(grating = 'GI', pointYaxis = py)
+    e4 = sl.caculateStarLightFilter(filter = 'i',pointYaxis = py)
+
+    e5=sl.caculateStrayLightFilter(filter = 'i', pixel_size_phy = 10, pixel_scale = 0.074, sun_pos = surveylist[i,6:9])
+    e6,_=sl.caculateStrayLightGrating(grating = 'GI', normFilter_es = 'g', sun_pos = surveylist[i,6:9])
 
     sky_pix[i,0] = e1
-    sky_pix[i,1] = e3
-    sky_pix[i,2] = e4
-    sky_pix[i,3] = v_mag
-    print(e1,e3,e4)
+    sky_pix[i,1] = e2
+    sky_pix[i,2] = e3
+    sky_pix[i,3] = e4
+    print(e1+e2,e1_+e3+e4,e5,e6)
 
     # print(e1,e2,e3,e4)