Commit e81eb82b authored by Yan Zhaojun's avatar Yan Zhaojun
Browse files

debug

parent ee15cc5b
Loading
Loading
Loading
Loading
Loading
+48 −48
Original line number Diff line number Diff line
@@ -2290,69 +2290,69 @@ class MCIsimulator():

# ########################################################################################################################################################################################################################################################

#     def zodiacal(self, ra, dec, time):
#         """
#             For given RA, DEC and TIME, return the interpolated zodical spectrum in Leinert-1998.
    def zodiacal(self, ra, dec, time):
        """
            For given RA, DEC and TIME, return the interpolated zodical spectrum in Leinert-1998.
    
#         :param ra: RA in unit of degree, ICRS frame
#         :param dec: DEC in unit of degree, ICRS frame
#         :param time: the specified string that in ISO format i.e., yyyy-mm-dd.
#         :return:
#             wave_A: wavelength of the zodical spectrum
#             spec_mjy: flux of the zodical spectrum, in unit of MJy/sr
#             spec_erg: flux of the zodical spectrum, in unit of erg/s/cm^2/A/sr
        :param ra: RA in unit of degree, ICRS frame
        :param dec: DEC in unit of degree, ICRS frame
        :param time: the specified string that in ISO format i.e., yyyy-mm-dd.
        :return:
            wave_A: wavelength of the zodical spectrum
            spec_mjy: flux of the zodical spectrum, in unit of MJy/sr
            spec_erg: flux of the zodical spectrum, in unit of erg/s/cm^2/A/sr
    
#         """
        """
    
#         # get solar position
#         dt = datetime.fromisoformat(time)
#         ###jd = julian.to_jd(dt, fmt='jd')
        # get solar position
        dt = datetime.fromisoformat(time)
        ###jd = julian.to_jd(dt, fmt='jd')
        
#         jd = time2jd(dt)
#         t = Time(jd, format='jd', scale='utc')
        jd = time2jd(dt)
        t = Time(jd, format='jd', scale='utc')
    
#         astro_sun = get_sun(t)
#         ra_sun, dec_sun = astro_sun.gcrs.ra.deg, astro_sun.gcrs.dec.deg
        astro_sun = get_sun(t)
        ra_sun, dec_sun = astro_sun.gcrs.ra.deg, astro_sun.gcrs.dec.deg
    
#         radec_sun = SkyCoord(ra=ra_sun*u.degree, dec=dec_sun*u.degree, frame='gcrs')
#         lb_sun = radec_sun.transform_to('geocentrictrueecliptic')
        radec_sun = SkyCoord(ra=ra_sun*u.degree, dec=dec_sun*u.degree, frame='gcrs')
        lb_sun = radec_sun.transform_to('geocentrictrueecliptic')
    
#         # get offsets between the target and sun.
#         radec_obj = SkyCoord(ra=ra*u.degree, dec=dec*u.degree, frame='icrs')
#         lb_obj = radec_obj.transform_to('geocentrictrueecliptic')
        # get offsets between the target and sun.
        radec_obj = SkyCoord(ra=ra*u.degree, dec=dec*u.degree, frame='icrs')
        lb_obj = radec_obj.transform_to('geocentrictrueecliptic')
    
#         beta = abs(lb_obj.lat.degree)
#         lamda = abs(lb_obj.lon.degree - lb_sun.lon.degree)
        beta = abs(lb_obj.lat.degree)
        lamda = abs(lb_obj.lon.degree - lb_sun.lon.degree)
    
#         # interpolated zodical surface brightness at 0.5 um
#         zodi = pd.read_csv(self.information['dir_path']+'MCI_inputData/refs/zodi_map.dat', sep='\s+', header=None, comment='#')
#         beta_angle = np.array([0, 5, 10, 15, 20, 25, 30, 45, 60, 75])
#         lamda_angle = np.array([0, 5, 10, 15, 20, 25, 30, 35, 40, 45,
#                               60, 75, 90, 105, 120, 135, 150, 165, 180])
#         xx, yy = np.meshgrid(beta_angle, lamda_angle)
#         f = interpolate.interp2d(xx, yy, zodi, kind='linear')
#         zodi_obj = f(beta, lamda)       # 10^�? W m�? sr�? um�?
        # interpolated zodical surface brightness at 0.5 um
        zodi = pd.read_csv(self.information['dir_path']+'MCI_inputData/refs/zodi_map.dat', sep='\s+', header=None, comment='#')
        beta_angle = np.array([0, 5, 10, 15, 20, 25, 30, 45, 60, 75])
        lamda_angle = np.array([0, 5, 10, 15, 20, 25, 30, 35, 40, 45,
                              60, 75, 90, 105, 120, 135, 150, 165, 180])
        xx, yy = np.meshgrid(beta_angle, lamda_angle)
        f = interpolate.interp2d(xx, yy, zodi, kind='linear')
        zodi_obj = f(beta, lamda)       # 10^�? W m�? sr�? um�?
    
#         # read the zodical spectrum in the ecliptic
#         cat_spec = pd.read_csv(self.information['dir_path']+'MCI_inputData/refs/solar_spec.dat', sep='\s+', header=None, comment='#')
#         wave = cat_spec[0].values       # A
#         spec0 = cat_spec[1].values      # 10^-8 W m^�? sr^�? μm^�?
#         zodi_norm = 252                 # 10^-8 W m^�? sr^�? μm^�?
        # read the zodical spectrum in the ecliptic
        cat_spec = pd.read_csv(self.information['dir_path']+'MCI_inputData/refs/solar_spec.dat', sep='\s+', header=None, comment='#')
        wave = cat_spec[0].values       # A
        spec0 = cat_spec[1].values      # 10^-8 W m^�? sr^�? μm^�?
        zodi_norm = 252                 # 10^-8 W m^�? sr^�? μm^�?
    
#         spec = spec0 * (zodi_obj / zodi_norm) * 1e-8  # W m^�? sr^�? μm^�?
        spec = spec0 * (zodi_obj / zodi_norm) * 1e-8  # W m^�? sr^�? μm^�?
    
#         # convert to the commonly used unit of MJy/sr, erg/s/cm^2/A/sr
#         wave_A = wave                                               # A
#         #spec_mjy = spec * 0.1 * wave_A**2 / 3e18 * 1e23 * 1e-6      # MJy/sr
#         spec_erg = spec * 0.1                                       # erg/s/cm^2/A/sr
#         spec_erg2 = spec_erg / 4.25452e10                           # erg/s/cm^2/A/arcsec^2
        # convert to the commonly used unit of MJy/sr, erg/s/cm^2/A/sr
        wave_A = wave                                               # A
        #spec_mjy = spec * 0.1 * wave_A**2 / 3e18 * 1e23 * 1e-6      # MJy/sr
        spec_erg = spec * 0.1                                       # erg/s/cm^2/A/sr
        spec_erg2 = spec_erg / 4.25452e10                           # erg/s/cm^2/A/arcsec^2
        
#         self.zodiacal_wave=wave_A    # in A
        self.zodiacal_wave=wave_A    # in A
        
#         self.zodiacal_flux=spec_erg2
        self.zodiacal_flux=spec_erg2
    
#         return wave_A, spec_erg2
    ###################################################################################
        return wave_A, spec_erg2
    ##################################################################################
    
    ##########################################################################
    def cal_Lensing_galaxy_img(self):