Upload New File

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

    # erg/s/cm^2/A/arcsec^2

    spec_erg2 = spec_erg / 4.25452e10

    return wave_A, spec_erg2

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

# path='/home/yan/MCI/'

# wave0, zodi0=zodiacal(10.0, 20.0, '2024-04-04', path)