Commit 4bb76b08 authored by Yan Zhaojun's avatar Yan Zhaojun
Browse files

Delete zodiacal.py

parent 5e1c7f2b
import numpy as np
import julian
from datetime import datetime
from astropy.time import Time
from astropy.coordinates import get_sun
from astropy.coordinates import SkyCoord
import pandas as pd
from astropy import units as u
from scipy import interpolate
def zodiacal(ra, dec, time, path):
"""
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.
:param path: the relative file path
:return:
wave_A: wavelength of the zodical spectrum
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')
###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
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')
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(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(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^�?
# 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
return wave_A, spec_erg2
###############################
# path='/home/yan/MCI/'
# wave0, zodi0=zodiacal(10.0, 20.0, '2024-04-04', path)
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment