Skip to content
GitLab
Commit a88bded7 authored by Zhang Xin's avatar Zhang Xin
Browse files

delete repeat file

parent a29e7df1
Hide whitespace changes
Inline Side-by-side
tests/TestSpecDisperse.py deleted 100644 → 0
View file @ a29e7df1
import unittest
from ObservationSim.MockObject.SpecDisperser import rotate90, SpecDisperser
from ObservationSim.Config import ChipOutput, Config
from ObservationSim.Instrument import Telescope, Chip, FilterParam, Filter, FocalPlane
from ObservationSim.MockObject import MockObject, Star
from ObservationSim.PSF import PSFGauss
import numpy as np
import galsim
from astropy.table import Table
from scipy import interpolate
import matplotlib.pyplot as plt
from lmfit.models import LinearModel, GaussianModel
from ObservationSim.Config.Header import generateExtensionHeader
import math
import yaml
def getAngle132(x1=0, y1=0, z1=0, x2=0, y2=0, z2=0, x3=0, y3=0, z3=0):
cosValue = 0;
angle = 0;
x11 = x1 - x3;
y11 = y1 - y3;
z11 = z1 - z3;
x22 = x2 - x3;
y22 = y2 - y3;
z22 = z2 - z3;
tt = np.sqrt((x11 * x11 + y11 * y11 + z11 * z11) * (x22 * x22 + y22 * y22 + z22 * z22));
if (tt == 0):
return 0;
cosValue = (x11 * x22 + y11 * y22 + z11 * z22) / tt;
if (cosValue > 1):
cosValue = 1;
if (cosValue < -1):
cosValue = -1;
angle = math.acos(cosValue);
return angle * 360 / (2 * math.pi);
def fit_SingleGauss(xX, yX, contmX, iHa0):
background = LinearModel(prefix='line_')
pars = background.make_params(intercept=yX.max(), slope=0)
pars = background.guess(yX, x=xX)
gauss = GaussianModel(prefix='g_')
pars.update(gauss.make_params())
pars['g_center'].set(iHa0, min=iHa0 - 3, max=iHa0 + 3)
pars['g_amplitude'].set(50, min=0)
pars['g_sigma'].set(12, min=0.0001)
mod = gauss + background
init = mod.eval(pars, x=xX)
outX = mod.fit(yX, pars, x=xX)
compsX = outX.eval_components(x=xX)
# print(outX.fit_report(min_correl=0.25))
# print outX.params['g_center']
outX.fit_report(min_correl=0.25)
# print(outX.fit_report(min_correl=0.25))
line_slopeX = float(outX.fit_report(min_correl=0.25).split('line_slope:')[1].split('+/-')[0]) * contmX
err_line_slopeX = float(
outX.fit_report(min_correl=0.25).split('line_slope:')[1].split('+/-')[1].split('(')[0]) * contmX
line_interceptX = float(outX.fit_report(min_correl=0.25).split('line_intercept:')[1].split('+/-')[0]) * contmX
err_line_interceptX = float(
outX.fit_report(min_correl=0.25).split('line_intercept:')[1].split('+/-')[1].split('(')[0]) * contmX
sigmaX = float(outX.fit_report(min_correl=0.25).split('g_sigma:')[1].split('+/-')[0])
err_sigmaX = float(outX.fit_report(min_correl=0.25).split('g_sigma:')[1].split('+/-')[1].split('(')[0])
fwhmX = float(outX.fit_report(min_correl=0.25).split('g_fwhm:')[1].split('+/-')[0])
err_fwhmX = float(outX.fit_report(min_correl=0.25).split('g_fwhm:')[1].split('+/-')[1].split('(')[0])
centerX = float(outX.fit_report(min_correl=0.25).split('g_center:')[1].split('+/-')[0])
err_centerX = float(outX.fit_report(min_correl=0.25).split('g_center:')[1].split('+/-')[1].split('(')[0])
return sigmaX, err_sigmaX, fwhmX, err_fwhmX, centerX, err_centerX
def produceObj(x,y,chip, ra, dec, pa):
pos_img = galsim.PositionD(chip.bound.xmin + x, chip.bound.ymin + y)
param = {}
param["star"] = 1
param["id"] = 1
param["ra"] = chip.img.wcs.posToWorld(pos_img).ra.deg
param["dec"] = chip.img.wcs.posToWorld(pos_img).dec.deg
param["z"] = 0
param["sed_type"] = 1
param["model_tag"] = 1
param["mag_use_normal"] = 10
obj = Star(param)
header_wcs = generateExtensionHeader(
xlen=chip.npix_x,
ylen=chip.npix_y,
ra=ra,
dec=dec,
pa=pa,
gain=chip.gain,
readout=chip.read_noise,
dark=chip.dark_noise,
saturation=90000,
psize=chip.pix_scale,
row_num=chip.rowID,
col_num=chip.colID,
extName='raw')
pos_img, offset, local_wcs, _ = obj.getPosImg_Offset_WCS(img=chip.img, chip=chip, img_header=header_wcs)
wave = np.arange(2500, 11000.5, 0.5)
# sedLen = wave.shape[0]
flux = pow(wave, -2) * 1e8
flux[200] = flux[200] * 10
flux[800] = flux[800] * 30
flux[2000] = flux[2000] * 5
obj.sed = Table(np.array([wave, flux]).T,
names=('WAVELENGTH', 'FLUX'))
return obj, pos_img
class TestSpecDisperse(unittest.TestCase):
def __init__(self, methodName='runTest',conff = '', throughputf = ''):
super(TestSpecDisperse,self).__init__(methodName)
self.conff = conff
self.throughputf = throughputf
def test_rotate901(self):
m = np.array([[1,2,3,4,5],[6,7,8,9,10],[11,12,13,14,15],[16,17,18,19,20],[21,22,23,24,25]])
m1 = np.array([[21,16,11,6,1],[22,17,12,7,2],[23,18,13,8,3],[24,19,14,9,4],[25,20,15,10,5]])
m2 = np.array([[5,10,15,20,25],[4,9,14,19,24],[3,8,13,18,23],[2,7,12,17,22],[1,6,11,16,21]])
xc = 2
yc = 2
isClockwise = 0
m1, xc1, yc1 = rotate90(array_orig=m, xc=xc, yc=yc, isClockwise=isClockwise)
self.assertTrue(xc1-xc == 0)
self.assertTrue(yc1-yc == 0)
self.assertTrue(np.sum(m-m1) == 0)
def test_rotate902(self):
m = np.array([[1,2,3,4,5],[6,7,8,9,10],[11,12,13,14,15],[16,17,18,19,20],[21,22,23,24,25]])
m1 = np.array([[21,16,11,6,1],[22,17,12,7,2],[23,18,13,8,3],[24,19,14,9,4],[25,20,15,10,5]])
m2 = np.array([[5,10,15,20,25],[4,9,14,19,24],[3,8,13,18,23],[2,7,12,17,22],[1,6,11,16,21]])
xc = 2
yc = 2
isClockwise =1
m1, xc1, yc1 = rotate90(array_orig=m, xc=xc, yc=yc, isClockwise=isClockwise)
self.assertTrue(xc1-xc == 0)
self.assertTrue(yc1-yc == 0)
self.assertTrue(np.sum(m-m2) == 0)
def test_Specdistperse1(self):
star = galsim.Gaussian(fwhm=0.39)
g_img = galsim.Image(100, 100, scale=0.074)
starImg = star.drawImage(image=g_img)
wave = np.arange(6200, 10000.5, 0.5)
# sedLen = wave.shape[0]
flux = pow(wave, -2) * 1e8
# flux[200] = flux[200] * 10
# flux[800] = flux[800] * 30
# flux[2000] = flux[2000] * 5
sed = Table(np.array([wave, flux]).T,
names=('WAVELENGTH', 'FLUX'))
conff = self.conff
throughput_f = self.throughputf
thp = Table.read(throughput_f)
thp_i = interpolate.interp1d(thp['WAVELENGTH'], thp['SENSITIVITY'])
sdp = SpecDisperser(orig_img=starImg, xcenter=0, ycenter=0, origin=[100, 100], tar_spec=sed, band_start=6200,
band_end=10000, isAlongY=0, conf=conff, gid=0)
spec = sdp.compute_spec_orders()
Aimg = spec['A'][0]
wave_pix = spec['A'][5]
wave_pos = spec['A'][3]
sens = spec['A'][6]
sh = Aimg.shape
spec_pix = np.zeros(sh[1])
for i in range(sh[1]):
spec_pix[i] = sum(Aimg[:, i])
# figure()
# imshow(Aimg)
wave_flux = np.zeros(wave_pix.shape[0])
for i in np.arange(1, wave_pix.shape[0] - 1):
w = wave_pix[i]
thp_w = thp_i(w)
deltW = (w - wave_pix[i - 1]) / 2 + (wave_pix[i + 1] - w) / 2
f = spec_pix[wave_pos[0] - 1 + i]
if 6200 <= w <= 10000:
f = f / thp_w
else:
f = 0
wave_flux[i] = f / deltW
sed_i = interpolate.interp1d(sed['WAVELENGTH'], sed['FLUX'])
ids = wave_pix < 9700
ids1 = wave_pix[ids] > 6500
print('Spec disperse flux test')
self.assertTrue(np.mean((wave_flux[ids][ids1] - sed_i(wave_pix[ids][ids1]))/sed_i(wave_pix[ids][ids1]))<0.004)
plt.figure()
plt.plot(wave_pix, wave_flux)
plt.plot(sed['WAVELENGTH'], sed['FLUX'])
plt.xlim(6200, 10000)
plt.ylim(1, 3)
plt.yscale('log')
plt.xlabel('$\lambda$')
plt.ylabel('$F\lambda$')
plt.legend(['extracted', 'input'])
plt.show()
def test_Specdistperse2(self):
psf_fwhm = 0.39
pix_scale = 0.074
star = galsim.Gaussian(fwhm=psf_fwhm)
g_img = galsim.Image(100, 100, scale=pix_scale)
starImg = star.drawImage(image=g_img)
wave = np.arange(6200, 10000.5, 0.5)
# sedLen = wave.shape[0]
flux = pow(wave, -2) * 1e8
flux[200] = flux[200] * 10
flux[800] = flux[800] * 30
flux[2000] = flux[2000] * 5
sed = Table(np.array([wave, flux]).T,
names=('WAVELENGTH', 'FLUX'))
conff = self.conff
throughput_f = self.throughputf
thp = Table.read(throughput_f)
thp_i = interpolate.interp1d(thp['WAVELENGTH'], thp['SENSITIVITY'])
sdp = SpecDisperser(orig_img=starImg, xcenter=0, ycenter=0, origin=[100, 100], tar_spec=sed, band_start=6200,
band_end=10000, isAlongY=0, conf=conff, gid=0)
spec = sdp.compute_spec_orders()
Aimg = spec['A'][0]
wave_pix = spec['A'][5]
wave_pos = spec['A'][3]
sens = spec['A'][6]
sh = Aimg.shape
spec_pix = np.zeros(sh[1])
for i in range(sh[1]):
spec_pix[i] = sum(Aimg[:, i])
# figure()
# imshow(Aimg)
wave_flux = np.zeros(wave_pix.shape[0])
for i in np.arange(1, wave_pix.shape[0] - 1):
w = wave_pix[i]
thp_w = thp_i(w)
deltW = (w - wave_pix[i - 1]) / 2 + (wave_pix[i + 1] - w) / 2
f = spec_pix[wave_pos[0] - 1 + i]
if 6200 <= w <= 10000:
f = f / thp_w
else:
f = 0
wave_flux[i] = f / deltW
sed_i = interpolate.interp1d(sed['WAVELENGTH'], sed['FLUX'])
input_em_lam = 6600
ids = wave_pix < input_em_lam+200
ids1 = wave_pix[ids] > input_em_lam-200
deltLamda_pix = (max(wave_pix[ids][ids1]) - min(wave_pix[ids][ids1])) / (wave_pix[ids][ids1].shape[0] - 1)
_, _, fwhmx, fwhmx_err, center, center_err = fit_SingleGauss(wave_pix[ids][ids1], wave_flux[ids][ids1], 1.0, 6600)
print('Emission line position and shape test')
self.assertTrue(input_em_lam-center < deltLamda_pix)
self.assertTrue(fwhmx/deltLamda_pix*pix_scale - psf_fwhm < np.abs(0.01))
# print('error is ',np.mean((wave_flux[ids][ids1] - sed_i(wave_pix[ids][ids1]))/sed_i(wave_pix[ids][ids1])))
# self.assertTrue(np.mean((wave_flux[ids][ids1] - sed_i(wave_pix[ids][ids1]))/sed_i(wave_pix[ids][ids1]))<0.004)
plt.figure()
plt.plot(wave_pix, wave_flux)
plt.plot(sed['WAVELENGTH'], sed['FLUX'])
plt.xlim(6200, 10000)
plt.ylim(1, 75)
plt.yscale('log')
plt.xlabel('$\lambda$')
plt.ylabel('$F\lambda$')
plt.legend(['extracted', 'input'])
plt.show()
def test_Specdistperse3(self):
psf_fwhm = 0.39
pix_scale = 0.074
star = galsim.Gaussian(fwhm=psf_fwhm)
g_img = galsim.Image(100, 100, scale=pix_scale)
starImg = star.drawImage(image=g_img)
wave = np.arange(6200, 10000.5, 0.5)
# sedLen = wave.shape[0]
flux = pow(wave, -2) * 1e8
flux[200] = flux[200] * 10
flux[800] = flux[800] * 30
flux[2000] = flux[2000] * 5
sed = Table(np.array([wave, flux]).T,
names=('WAVELENGTH', 'FLUX'))
conff = self.conff
throughput_f = self.throughputf
thp = Table.read(throughput_f)
thp_i = interpolate.interp1d(thp['WAVELENGTH'], thp['SENSITIVITY'])
sdp = SpecDisperser(orig_img=starImg, xcenter=0, ycenter=0, origin=[100, 100], tar_spec=sed, band_start=6200,
band_end=8000, isAlongY=0, conf=conff, gid=0)
sdp1 = SpecDisperser(orig_img=starImg, xcenter=0, ycenter=0, origin=[100, 100], tar_spec=sed, band_start=8000,
band_end=10000, isAlongY=0, conf=conff, gid=0)
spec = sdp.compute_spec_orders()
spec1 = sdp1.compute_spec_orders()
Aimg = spec['A'][0] + spec1['A'][0]
wave_pix = spec['A'][5]
wave_pos = spec['A'][3]
sens = spec['A'][6]
sh = Aimg.shape
spec_pix = np.zeros(sh[1])
for i in range(sh[1]):
spec_pix[i] = sum(Aimg[:, i])
wave_flux = np.zeros(wave_pix.shape[0])
for i in np.arange(1, wave_pix.shape[0] - 1):
w = wave_pix[i]
thp_w = thp_i(w)
deltW = (w - wave_pix[i - 1]) / 2 + (wave_pix[i + 1] - w) / 2
f = spec_pix[wave_pos[0] - 1 + i]
if 6200 <= w <= 10000:
f = f / thp_w
else:
f = 0
wave_flux[i] = f / deltW
sdp2 = SpecDisperser(orig_img=starImg, xcenter=0, ycenter=0, origin=[100, 100], tar_spec=sed, band_start=6200,
band_end=10000, isAlongY=0, conf=conff, gid=0)
spec2 = sdp2.compute_spec_orders()
Aimg2 = spec2['A'][0]
spec_pix2 = np.zeros(sh[1])
for i in range(sh[1]):
spec_pix2[i] = sum(Aimg2[:, i])
wave_flux2 = np.zeros(wave_pix.shape[0])
for i in np.arange(1, wave_pix.shape[0] - 1):
w = wave_pix[i]
thp_w = thp_i(w)
deltW = (w - wave_pix[i - 1]) / 2 + (wave_pix[i + 1] - w) / 2
f = spec_pix2[wave_pos[0] - 1 + i]
if 6200 <= w <= 10000:
f = f / thp_w
else:
f = 0
wave_flux2[i] = f / deltW
r1_i = interpolate.interp1d(wave_pix, wave_flux2)
r2_i = interpolate.interp1d(wave_pix, wave_flux)
print('Spec Splicing test')
self.assertTrue(r1_i(8000)-r2_i(8000) < np.abs(0.0001))
# self.assertTrue(fwhmx/deltLamda_pix*pix_scale - psf_fwhm < np.abs(0.01))
# print('error is ',np.mean((wave_flux[ids][ids1] - sed_i(wave_pix[ids][ids1]))/sed_i(wave_pix[ids][ids1])))
# self.assertTrue(np.mean((wave_flux[ids][ids1] - sed_i(wave_pix[ids][ids1]))/sed_i(wave_pix[ids][ids1]))<0.004)
plt.figure()
plt.plot(wave_pix, wave_flux2)
plt.plot(wave_pix, wave_flux)
# plt.plot(sed['WAVELENGTH'], sed['FLUX'])
plt.xlim(6200, 10000)
plt.ylim(1, 4)
plt.yscale('log')
plt.xlabel('$\lambda$')
plt.ylabel('$F\lambda$')
plt.legend(['one spec', 'split in 8000 A'])
plt.show()
def test_double_disperse(self):
work_dir = "/public/home/fangyuedong/CSST_unittest/CSST/test/"
# data_dir = "/Volumes/Extreme SSD/SimData/"
data_dir = "/data/simudata/CSSOSDataProductsSims/data/"
configFn = '/Users/zhangxin/Work/SlitlessSim/CSST_SIM/CSST_C6/csst-simulation/config/config_C6.yaml'
normFilterFn = "/Users/zhangxin/Work/SlitlessSim/CSST_SIM/CSST_C6/csst-simulation/Catalog/data/SLOAN_SDSS.g.fits"
norm_star = Table.read(normFilterFn)
with open(configFn, "r") as stream:
try:
config = yaml.safe_load(stream)
for key, value in config.items():
print(key + " : " + str(value))
except yaml.YAMLError as exc:
print(exc)
# config = Config.read_config(configFn)
# path_dict = Config.config_dir(config,work_dir=work_dir, data_dir=data_dir)
filter_param = FilterParam()
focal_plane = FocalPlane(survey_type=config["obs_setting"]["survey_type"])
chip = Chip(1, config=config)
filter_id, filter_type = chip.getChipFilter()
filt = Filter(filter_id=filter_id, filter_type=filter_type, filter_param=filter_param,
ccd_bandpass=chip.effCurve)
tel = Telescope()
psf_model = PSFGauss(chip=chip)
wcs_fp = focal_plane.getTanWCS(float(config["obs_setting"]["ra_center"]), float(config["obs_setting"]["dec_center"]), float(config["obs_setting"]["image_rot"]) * galsim.degrees, chip.pix_scale)
chip.img = galsim.ImageF(chip.npix_x, chip.npix_y)
chip.img.setOrigin(chip.bound.xmin, chip.bound.ymin)
chip.img.wcs = wcs_fp
obj, pos_img = produceObj(2000,4500, chip,float(config["obs_setting"]["ra_center"]), float(config["obs_setting"]["dec_center"]), float(config["obs_setting"]["image_rot"]))
obj.drawObj_slitless(
tel=tel,
pos_img=pos_img,
psf_model=psf_model,
bandpass_list=filt.bandpass_sub_list,
filt=filt,
chip=chip,
g1=0,
g2=0,
exptime=150,
normFilter=norm_star)
obj, pos_img = produceObj(3685, 6500, chip,float(config["obs_setting"]["ra_center"]), float(config["obs_setting"]["dec_center"]), float(config["obs_setting"]["image_rot"]))
obj.drawObj_slitless(
tel=tel,
pos_img=pos_img,
psf_model=psf_model,
bandpass_list=filt.bandpass_sub_list,
filt=filt,
chip=chip,
g1=0,
g2=0,
exptime=150,
normFilter=norm_star)
obj, pos_img = produceObj(5000, 2500, chip, float(config["obs_setting"]["ra_center"]), float(config["obs_setting"]["dec_center"]), float(config["obs_setting"]["image_rot"]))
obj.drawObj_slitless(
tel=tel,
pos_img=pos_img,
psf_model=psf_model,
bandpass_list=filt.bandpass_sub_list,
filt=filt,
chip=chip,
g1=0,
g2=0,
exptime=150,
normFilter=norm_star)
print('Spec double disperse test')
from astropy.io import fits
fits.writeto('test.fits',chip.img.array, overwrite = True)
# plt.figure()
# plt.imshow(chip.img.array)
# plt.show()
def test_SLSImage_rotation(self):
from astropy.wcs import WCS
configFn = '/Users/zhangxin/Work/SlitlessSim/CSST_SIM/CSST_C6/csst-simulation/config/config_C6.yaml'
with open(configFn, "r") as stream:
try:
config = yaml.safe_load(stream)
for key, value in config.items():
print(key + " : " + str(value))
except yaml.YAMLError as exc:
print(exc)
chip = Chip(1, config=config)
ra=float(config["obs_setting"]["ra_center"])
dec=float(config["obs_setting"]["dec_center"])
pa=float(config["obs_setting"]["image_rot"])
header_wcs1 = generateExtensionHeader(
xlen=chip.npix_x,
ylen=chip.npix_y,
ra=ra,
dec=dec,
pa=pa,
gain=chip.gain,
readout=chip.read_noise,
dark=chip.dark_noise,
saturation=90000,
psize=chip.pix_scale,
row_num=chip.rowID,
col_num=chip.colID,
extName='raw', center_rot=0)
center = np.array([chip.npix_x / 2, chip.npix_y / 2])
h_wcs1 = WCS(header_wcs1)
x1, y1 = center + [100,0]
sky_1 = h_wcs1.pixel_to_world(x1,y1)
rot_angle = 1
header_wcs2 = generateExtensionHeader(
xlen=chip.npix_x,
ylen=chip.npix_y,
ra=ra,
dec=dec,
pa=pa,
gain=chip.gain,
readout=chip.read_noise,
dark=chip.dark_noise,
saturation=90000,
psize=chip.pix_scale,
row_num=chip.rowID,
col_num=chip.colID,
extName='raw', center_rot=rot_angle)
h_wcs2 = WCS(header_wcs2)
x2, y2 = h_wcs2.world_to_pixel(sky_1)
angle = getAngle132(x1,y1,0,x2,y2,0,center[0],center[1],0)
print(angle)
self.assertTrue(rot_angle - angle < np.abs(0.001))
rot_angle = 10
header_wcs2 = generateExtensionHeader(
xlen=chip.npix_x,
ylen=chip.npix_y,
ra=ra,
dec=dec,
pa=pa,
gain=chip.gain,
readout=chip.read_noise,
dark=chip.dark_noise,
saturation=90000,
psize=chip.pix_scale,
row_num=chip.rowID,
col_num=chip.colID,
extName='raw', center_rot=rot_angle)
h_wcs2 = WCS(header_wcs2)
x2, y2 = h_wcs2.world_to_pixel(sky_1)
angle = getAngle132(x1, y1, 0, x2, y2, 0, center[0], center[1], 0)
print(angle)
self.assertTrue(rot_angle - angle < np.abs(0.001))
rot_angle = 50
header_wcs2 = generateExtensionHeader(
xlen=chip.npix_x,
ylen=chip.npix_y,
ra=ra,
dec=dec,
pa=pa,
gain=chip.gain,
readout=chip.read_noise,
dark=chip.dark_noise,
saturation=90000,
psize=chip.pix_scale,
row_num=chip.rowID,
col_num=chip.colID,
extName='raw', center_rot=rot_angle)
h_wcs2 = WCS(header_wcs2)
x2, y2 = h_wcs2.world_to_pixel(sky_1)
angle = getAngle132(x1, y1, 0, x2, y2, 0, center[0], center[1], 0)
print(angle)
self.assertTrue(rot_angle - angle < np.abs(0.001))
chip = Chip(27, config=config)
ra = float(config["obs_setting"]["ra_center"])
dec = float(config["obs_setting"]["dec_center"])
pa = float(config["obs_setting"]["image_rot"])
header_wcs1 = generateExtensionHeader(
xlen=chip.npix_x,
ylen=chip.npix_y,
ra=ra,
dec=dec,
pa=pa,
gain=chip.gain,
readout=chip.read_noise,
dark=chip.dark_noise,
saturation=90000,
psize=chip.pix_scale,
row_num=chip.rowID,
col_num=chip.colID,
extName='raw', center_rot=0)
center = np.array([chip.npix_x / 2, chip.npix_y / 2])
h_wcs1 = WCS(header_wcs1)
x1, y1 = center + [100, 0]
sky_1 = h_wcs1.pixel_to_world(x1, y1)
rot_angle = 1
header_wcs2 = generateExtensionHeader(
xlen=chip.npix_x,
ylen=chip.npix_y,
ra=ra,
dec=dec,
pa=pa,
gain=chip.gain,
readout=chip.read_noise,
dark=chip.dark_noise,
saturation=90000,
psize=chip.pix_scale,
row_num=chip.rowID,
col_num=chip.colID,
extName='raw', center_rot=rot_angle)
h_wcs2 = WCS(header_wcs2)
x2, y2 = h_wcs2.world_to_pixel(sky_1)
angle = getAngle132(x1, y1, 0, x2, y2, 0, center[0], center[1], 0)
print(angle)
self.assertTrue(rot_angle - angle < np.abs(0.001))
rot_angle = 10
header_wcs2 = generateExtensionHeader(
xlen=chip.npix_x,
ylen=chip.npix_y,
ra=ra,
dec=dec,
pa=pa,
gain=chip.gain,
readout=chip.read_noise,
dark=chip.dark_noise,
saturation=90000,
psize=chip.pix_scale,
row_num=chip.rowID,
col_num=chip.colID,
extName='raw', center_rot=rot_angle)
h_wcs2 = WCS(header_wcs2)
x2, y2 = h_wcs2.world_to_pixel(sky_1)
angle = getAngle132(x1, y1, 0, x2, y2, 0, center[0], center[1], 0)
print(angle)
self.assertTrue(rot_angle - angle < np.abs(0.001))
rot_angle = 50
header_wcs2 = generateExtensionHeader(
xlen=chip.npix_x,
ylen=chip.npix_y,
ra=ra,
dec=dec,
pa=pa,
gain=chip.gain,
readout=chip.read_noise,
dark=chip.dark_noise,
saturation=90000,
psize=chip.pix_scale,
row_num=chip.rowID,
col_num=chip.colID,
extName='raw', center_rot=rot_angle)
h_wcs2 = WCS(header_wcs2)
x2, y2 = h_wcs2.world_to_pixel(sky_1)
angle = getAngle132(x1, y1, 0, x2, y2, 0, center[0], center[1], 0)
print(angle)
self.assertTrue(rot_angle - angle < np.abs(0.001))
if __name__ == '__main__':
conff='/Users/zhangxin/Work/SlitlessSim/CSST_SIM/CSST_C6/csst-simulation/ObservationSim/Instrument/data/sls_conf/CSST_GI2.conf'
throughputf='/Users/zhangxin/Work/SlitlessSim/CSST_SIM/CSST_C6/csst-simulation/ObservationSim/Instrument/data/sls_conf/GI.Throughput.1st.fits'
suit = unittest.TestSuite()
case1 = TestSpecDisperse('test_Specdistperse1',conff,throughputf)
suit.addTest(case1)
case2 = TestSpecDisperse('test_Specdistperse2', conff, throughputf)
suit.addTest(case2)
case3 = TestSpecDisperse('test_Specdistperse3', conff, throughputf)
suit.addTest(case3)
case4 = TestSpecDisperse('test_double_disperse', conff, throughputf)
suit.addTest(case4)
case5 = TestSpecDisperse('test_SLSImage_rotation')
suit.addTest(case5)
unittest.TextTestRunner(verbosity=2).run(suit)
# runner = unittest.TextTestRunner()
# runner.run(suit)
\ No newline at end of file
tests/TestStraylight.py deleted 100644 → 0
View file @ a29e7df1
import unittest
from ObservationSim.Straylight import Straylight
import numpy as np
import math
import astropy.constants as cons
import galsim
from astropy.table import Table
from scipy import interpolate
import matplotlib.pyplot as plt
hubbleAverZodiacal = {'nuv':0.0035,'u':0.0163,'g':0.1109,'r':0.1471,'i':0.1568,'z':0.0953,'y':0.0283}
hubbleAverEarthShine = {'nuv':0.00024,'u':0.0051,'g':0.0506,'r':0.0591,'i':0.0568,'z':0.0315,'y':0.0090}
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);
z1 = np.sin(dec / 57.2957795);
return np.array([x1, y1, z1])
def getAngle132(x1=0, y1=0, z1=0, x2=0, y2=0, z2=0, x3=0, y3=0, z3=0):
cosValue = 0;
angle = 0;
x11 = x1 - x3;
y11 = y1 - y3;
z11 = z1 - z3;
x22 = x2 - x3;
y22 = y2 - y3;
z22 = z2 - z3;
tt = np.sqrt((x11 * x11 + y11 * y11 + z11 * z11) * (x22 * x22 + y22 * y22 + z22 * z22));
if (tt == 0):
return 0;
cosValue = (x11 * x22 + y11 * y22 + z11 * z22) / tt;
if (cosValue > 1):
cosValue = 1;
if (cosValue < -1):
cosValue = -1;
angle = math.acos(cosValue);
return angle * 360 / (2 * math.pi);
def calculateAnglePwithEarth(sat = np.array([0,0,0]), pointing = np.array([0,0,0]), sun = np.array([0,0,0])):
modSat = np.sqrt(sat[0]*sat[0] + sat[1]*sat[1]+sat[2]*sat[2])
modPoint = np.sqrt(pointing[0]*pointing[0] + pointing[1]*pointing[1] + pointing[2]*pointing[2])
withLocalZenithAngle = (pointing[0] * sat[0] + pointing[1] * sat[1] + pointing[2] * sat[2]) / (modPoint*modSat)
innerM_sat_sun = sat[0] * sun[0] + sat[1] * sun[1] + sat[2] * sun[2]
cosAngle = innerM_sat_sun / (modSat * cons.au.value/1000)
isInSunSide = 1
if (cosAngle < -0.3385737): #cos109.79
isInSunSide = -1;
elif cosAngle >= -0.3385737 and cosAngle <= 0.3385737:
isInSunSide = 0;
return math.acos(withLocalZenithAngle)*180/math.pi,isInSunSide
class TestStraylight(unittest.TestCase):
def __init__(self, methodName='runTest',datFn = '', filter = 'i', grating = "GI"):
super(TestStraylight,self).__init__(methodName)
self.pointingData = np.loadtxt(datFn, dtype=np.double)
self.filter = filter
self.grating = grating
def test_EarthShineFilter(self):
d_sh = self.pointingData.shape
sl_e_pix = np.zeros([d_sh[0],3],dtype=np.double)
for i in np.arange(d_sh[0]):
# if i > 50:
# continue
ju = self.pointingData[i, 5]
# pointing = transRaDec2D(self.pointingData[i, 0], self.pointingData[i, 1])
# print(ju, pointing, surveylist[i,3:9])
sl = Straylight(jtime=ju, sat_pos=self.pointingData[i, 6:9], pointing_radec=np.array([self.pointingData[i, 0], self.pointingData[i, 1]]),sun_pos=self.pointingData[i,9:12])
e1, py = sl.calculateEarthShineFilter(filter=self.filter)
earthZenithAngle, isInSunSide = calculateAnglePwithEarth(sat=self.pointingData[i, 6:9], pointing= sl.pointing, sun=self.pointingData[i,9:12])
# e2, _ = sl.calculateZodiacalFilter2(filter='i', sun_pos=sl.sun_pos)
# e3 = sl.calculateStarLightFilter(filter='i', pointYaxis=py)
# e_all = sl.calculateStrayLightFilter(filter='i')
# s_pix, spec = sl.calculateStrayLightGrating(grating='GI')
sl_e_pix[i,0] = e1
sl_e_pix[i, 1] = earthZenithAngle
sl_e_pix[i, 2] = isInSunSide
median = np.median(sl_e_pix[:,0])
print(' average Earthshine %s: %e' % (self.filter, median))
self.assertTrue(median-hubbleAverEarthShine[self.filter] < 0.1)
plt.figure()
ids1 = sl_e_pix[:, 2] == 1
ids2 = sl_e_pix[:, 2] != 1
plt.plot(sl_e_pix[ids1, 0], sl_e_pix[ids1, 1], 'r.')
plt.plot(sl_e_pix[ids2, 0], sl_e_pix[ids2, 1], 'b.')
plt.legend(['In Sun Side', 'In Earths shadow'])
plt.xlabel('straylight-earthshine(e-/pixel/s)')
plt.ylabel('Angle with local zenith(degree)')
plt.show()
def test_ZodiacalFilter(self):
d_sh = self.pointingData.shape
sl_e_pix = np.zeros([d_sh[0],2],dtype=np.double)
for i in np.arange(d_sh[0]):
ju = self.pointingData[i, 5]
sl = Straylight(jtime=ju, sat_pos=self.pointingData[i, 6:9], pointing_radec=np.array([self.pointingData[i, 0], self.pointingData[i, 1]]),sun_pos=self.pointingData[i,9:12])
e1, _ = sl.calculateZodiacalFilter2(filter=self.filter, sun_pos=sl.sun_pos)
sl_e_pix[i,0] = e1
sl_e_pix[i,1] = getAngle132(x1=self.pointingData[i,9], y1=self.pointingData[i,10], z1=self.pointingData[i,11], x2=sl.pointing[0],
y2=sl.pointing[1], z2=sl.pointing[2], x3=0, y3=0, z3=0)
plt.figure()
plt.plot(sl_e_pix[:, 0], sl_e_pix[:, 1], 'r.')
plt.xlabel('straylight-zodiacal(e-/pixel/s)')
plt.ylabel('Angle between pointing and sun(degree)')
plt.show()
median = np.median(sl_e_pix[:,0])
print(' average Zodiacal %s: %f' % (self.filter, median))
self.assertTrue(median-hubbleAverZodiacal[self.filter] < 0.1)
def test_StarFilter(self):
d_sh = self.pointingData.shape
sl_e_pix = np.zeros(d_sh[0],dtype=np.double)
tnum = 10
for i in np.arange(tnum):
# if i > 50:
# continue
ju = self.pointingData[i, 5]
# pointing = transRaDec2D(self.pointingData[i, 0], self.pointingData[i, 1])
# print(ju, pointing, surveylist[i,3:9])
sl = Straylight(jtime=ju, sat_pos=self.pointingData[i, 6:9], pointing_radec=np.array([self.pointingData[i, 0], self.pointingData[i, 1]]),sun_pos=self.pointingData[i,9:12])
e1, py = sl.calculateEarthShineFilter(filter=self.filter)
# e2, _ = sl.calculateZodiacalFilter2(filter='i', sun_pos=sl.sun_pos)
e3 = sl.calculateStarLightFilter(filter=self.filter, pointYaxis=py)
# e_all = sl.calculateStrayLightFilter(filter='i')
# s_pix, spec = sl.calculateStrayLightGrating(grating='GI')
sl_e_pix[i] = e3
median = np.median(sl_e_pix[0:tnum])
print(' average Earthshine %s: %e' % (self.filter, median))
self.assertTrue(median-hubbleAverEarthShine[self.filter] < 0.2)
def test_GratingStraylight(self):
d_sh = self.pointingData.shape
sl_e_pix = np.zeros(d_sh[0],dtype=np.double)
tnum = 10
for i in np.arange(tnum):
# if i > 50:
# continue
ju = self.pointingData[i, 5]
# pointing = transRaDec2D(self.pointingData[i, 0], self.pointingData[i, 1])
# print(ju, pointing, surveylist[i,3:9])
sl = Straylight(jtime=ju, sat_pos=self.pointingData[i, 6:9], pointing_radec=np.array([self.pointingData[i, 0], self.pointingData[i, 1]]),sun_pos=self.pointingData[i,9:12])
# e1, py = sl.calculateEarthShineFilter(filter=self.filter)
# e2, _ = sl.calculateZodiacalFilter2(filter='i', sun_pos=sl.sun_pos)
# e3 = sl.calculateStarLightFilter(filter=self.filter, pointYaxis=py)
# e_all = sl.calculateStrayLightFilter(filter='i')
s_pix, spec = sl.calculateStrayLightGrating(grating=self.grating)
sl_e_pix[i] = s_pix
plt.figure()
plt.plot(spec['WAVELENGTH'], spec['FLUX'], 'r')
plt.xlabel('WAVELENGTH')
plt.ylabel('F$\lambda$(erg/s/cm2/A/arcsec2)')
plt.xlim(2000,10000)
plt.show()
median = np.median(sl_e_pix[0:tnum])
print(' average Earthshine %s: %e' % (self.grating, median))
self.assertTrue(median < 0.8)
if __name__ == '__main__':
suit = unittest.TestSuite()
# case1 = TestStraylight('test_EarthShineFilter',datFn = 'Straylight_test.dat', filter = 'i')
# suit.addTest(case1)
# case2 = TestStraylight('test_ZodiacalFilter',datFn = 'Straylight_test.dat',filter = 'i')
# suit.addTest(case2)
# case3 = TestStraylight('test_StarFilter', datFn='Straylight_test.dat', filter='i')
# suit.addTest(case3)
case4 = TestStraylight('test_GratingStraylight', datFn='Straylight_test.dat', grating = 'GI')
suit.addTest(case4)
unittest.TextTestRunner(verbosity=2).run(suit)
\ No newline at end of file
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