diff --git a/scopesim/effects/__init__.py b/scopesim/effects/__init__.py index e232667f..59658d69 100644 --- a/scopesim/effects/__init__.py +++ b/scopesim/effects/__init__.py @@ -17,5 +17,8 @@ from .fits_headers import * from .rotation import * +from .mosaic_atmo_disp_fibre_coupling import MosaicAtmosDispFibreCoupling +from . import mosaic_atmo_disp_fibre_coupling_utils + # from . import effects_utils diff --git a/scopesim/effects/mosaic_atmo_disp_fibre_coupling.py b/scopesim/effects/mosaic_atmo_disp_fibre_coupling.py new file mode 100644 index 00000000..744379ff --- /dev/null +++ b/scopesim/effects/mosaic_atmo_disp_fibre_coupling.py @@ -0,0 +1,321 @@ +import numpy as np +from astropy import units as u +from astropy.coordinates import Angle +import matplotlib.pyplot as plt +import matplotlib as mpl + +from . import mosaic_atmo_disp_fibre_coupling_utils as disp_utils +from . import TERCurve + + +class MosaicAtmosDispFibreCoupling(TERCurve): + def __init__(self, **kwargs): + params = {"HA_start": 0, + "HA_end": 0, + "declination": 0, + "band": "HR_NIR_H", + "sampling": 0.01, # um + "guide_waveref": -1, + "aperture_waveref": -1} + + super(TERCurve, self).__init__() + self.meta.update(params) + self.meta.update(kwargs) + + self.coupling_fractions = AtmosDispFibreCoupling(**kwargs) + + def apply_to(self, obj, **kwargs): + + if isinstance(obj, SourceBase): + # return objects are: 1D wavelength, 2D transmission curves + waves, trans_arr = self.coupling_fractions.run() + + + return base + + +class AtmosDispFibreCoupling: + """ + Class to quantify atmospheric dispersion effects on a MOSAIC integration. + See .run for inputs and outputs + + Author: Jay Stephan + + Bugs/Issues: + HR Wavelength and aperture major axis values need to be confirmed + """ + def __init__(self, **kwargs): + #Simulation parameters for MOSAIC + params = {'telescope_diameter':39, #m, diameter of telescope (ELT) + 'wavefront_outerscale':46, #m, wavefront outer scale for FWHM change with airmass/wavelength. Value to be confirmed. + 'median_seeing':.68, #arcsec, median seeing at Paranal + 'median_seeing_wl':.5, #um, wavelength median seeing corresponds to + + 'latitude':-24.6272, #deg, MUST be negative (southern hemisphere) for now + 'temperature':10+273.15, #K, temperature at Paranal + 'humidity':14.5/100, #fraction, relative humidity at Paranal + 'pressure':750, #mbar, pressure at Paranal + + 'LR_VIS_major_axis':.702, #arcsec, major axis of observing modes apertures + 'LR_VIS_hexagon_radius': 2, #radius of aperture hexagon array in hexagons + 'LR_NIR_major_axis':.57, + 'LR_NIR_hexagon_radius': 2, + 'HR_VIS_major_axis':.700, + 'HR_VIS_hexagon_radius': 3, + 'HR_NIR_major_axis':.57, #This needs to be confirmed, but should be correct + 'HR_NIR_hexagon_radius': 2, + + 'custom_hexagon_radius': 0, #Change to be non-zero to overide aperture hexagon array + + 'LR_VIS_B':[.390,.458], #um, MOSAIC bands + 'LR_VIS_G':[.450,.591], + 'LR_VIS_R':[.586,.770], + 'LR_VIS_All':[.390,.770], + 'LR_NIR_IY':[.770,1.063], + 'LR_NIR_J':[1.01,1.395], + 'LR_NIR_H':[1.420,1.857], + 'LR_NIR_All':[.770,1.857], + 'HR_VIS_G':[.510,.568], #HR values need to be confirmed, especially VIS + 'HR_VIS_R':[.610,.680], + 'HR_NIR_IY':[.770,.907], + 'HR_NIR_H':[1.523,1.658], + + 'sim_scale':.005, #arcsec/pixel, scale to carry out the simulation - smaller is slower, more accurate. Do not put above 0.01 + 'sim_HA_samples':21, #number of instantaneous snapshots to average over for the integration + 'relative_plate_PA_angle':0, #deg, relative angle of the plate/apertures and PA=0. For PA=0 along semi major axis, set to 90 deg + } + + # for kwarg in kwargs.items(): + # params[kwarg[0]]=kwarg[1] + params.update(kwargs) + self.config = params + self.input={} + self.output={} + + def load_MOSAIC_band(self,band,sampling): + """ + Generates wavelengths for the simulation from the chosen band, + and acquires relevant diameter for the band (changes between LR/HR, and VIS/NIR) + """ + self.input['sampling']=sampling + self.input['band']=band + + wave_min,wave_max= self.config[band][0],self.config[band][1] + + self.output['wavelengths'] = np.arange(wave_min,wave_max,sampling) + self.output['major_axis']=self.config[band[0:6]+"_major_axis"] + + def load_HA(self,HA_start,HA_end,declination): + """ + Calculates the airmasses to use in the simulation using the observation parameters (HA and dec) + """ + HA_range=np.linspace(HA_start,HA_end,self.config['sim_HA_samples']) + self.input['HA_range']=HA_range + self.input['targ_dec']=declination + + #latitude needs to be negative for now + lat = self.config['latitude'] * np.pi/180 + dec = declination*np.pi/180 + + #Need to check if the target is below the horizon for the given list of HA, and if so remove it. + LHA_below_horizon=np.rad2deg(np.arccos(-np.tan(lat)*np.tan(dec)))/15 + if str(LHA_below_horizon) != 'nan': + for val in HA_range.copy(): + if abs(val) > abs(LHA_below_horizon): + print("At HA %2.2fh, target goes below horizon - removing this from HA range" % (val)) + HA_range.remove(val) + if dec > np.pi/2 + lat: #If the target has a too high declination, it will never be seen at Cerro Paranal + print("Target always below Horizon") + return + + airmasses=1/(np.sin(lat)*np.sin(dec)+np.cos(lat)*np.cos(dec)*np.cos(Angle(HA_range*u.hour).rad)) + self.output['airmasses']=np.array(airmasses) + + para_angles=disp_utils.parallatic_angle(np.array(HA_range),dec,lat) + self.output['raw_para_angles']=np.array(para_angles) #actual PAs + + def calculate_shifts(self, guide_waveref, aperture_waveref): + """ + Calculates shifts of the wavelengths at the different airmasses using Fillipenko + """ + self.input['guide_waveref']=guide_waveref + self.input['aperture_waveref']=aperture_waveref + + airmasses=self.output['airmasses'] + wavelengths=self.output['wavelengths'] + + #centring refers to the index of the hour angles at which we centre the aperture/guiding on a wavelength + centring_index=int((len(airmasses)-1)/2) + + centre_shift=disp_utils.diff_shift(aperture_waveref,airmasses[centring_index],guide_waveref,self.config) #shift of the original aperture centre wavelength from guide wavelength + centring_q=self.output['raw_para_angles'][centring_index] + + raw_para_angles=self.output['raw_para_angles'] + para_angles=self.output['raw_para_angles'].copy() + for i in range(0,len(para_angles)): #change in PAs from centring index + para_angles[i]=para_angles[i]-self.output['raw_para_angles'][centring_index] + + shifts_para=[] + phi=np.deg2rad(self.config['relative_plate_PA_angle']) + for count,airmass in enumerate(airmasses): #for each airmass, calculate AD shift + shift_vals=disp_utils.diff_shift(wavelengths,airmass,guide_waveref,self.config) + airmass_shifts=[] + + for i in range(0,len(shift_vals)): + x=(shift_vals[i])*np.sin(raw_para_angles[count])-centre_shift*np.sin(centring_q) + y=(shift_vals[i])*np.cos(raw_para_angles[count])-centre_shift*np.cos(centring_q) + airmass_shifts.append([x*np.cos(phi)-y*np.sin(phi),y*np.cos(phi)+x*np.sin(phi)]) + + shifts_para.append(airmass_shifts) + + self.output['shifts']=np.array(shifts_para) + centre_shift_para=[-centre_shift*np.sin(centring_q),-centre_shift*np.cos(centring_q)] + centre_shift_para=[centre_shift_para[0]*np.cos(phi)-centre_shift_para[1]*np.sin(phi), + centre_shift_para[1]*np.cos(phi)+centre_shift_para[0]*np.sin(phi)] + self.output['centre_shift']=centre_shift_para + + def load_PSFs(self): + """ + Generates the shifted moffat PSFs - one for for wavelength at each airmass + """ + all_PSFs=[] + + for i in range(0,len(self.output['airmasses'])): + PSFs=disp_utils.make_moffat_PSFs(self.output['wavelengths'],self.output['shifts'][i], + self.output['airmasses'][i],self.output['major_axis'],self.config) + all_PSFs.append(PSFs) + + self.output['PSFs']=all_PSFs + + def load_aperture(self,band): + """ + Generates apertures, one for each fibre in the bundle + """ + apertures,apertures_table=disp_utils.make_aperture(band,self.output['major_axis'],self.config) + self.output['apertures']=apertures + self.apertures_table=apertures_table + + def calculate_transmissions(self): + """ + Calculates the transmission of the PSFs through the apertures. + """ + apertures_resized = np.repeat(self.output['apertures'][:,np.newaxis], np.shape(self.output['PSFs'])[1], axis=1) + apertures_resized = np.repeat(apertures_resized[:,np.newaxis], np.shape(self.output['PSFs'])[0], axis=1) + + PSFs_through_apertures=apertures_resized*self.output['PSFs'] + self.output['PSFs_through_apertures']=PSFs_through_apertures + + transmissions=np.sum(PSFs_through_apertures,axis=(-1,-2)) + self.output['raw_transmissions']=transmissions + + integration_transmissions=np.mean(transmissions,axis=1) + self.output['raw_integration_transmissions']=integration_transmissions #collapsed along airmass axis + + def plots(self): + """ + Function to illustrate simulation results + 1) Transmission vs wavelength curves for individual fibres and entire bundle + 2) Track plot of monochromatic spot PSFs on the aperture over an integration + """ + plt.style.use('bmh') + + fig,ax=plt.subplots(figsize=[7,5]) + plt.ylim(0,0.1) + ax.set_ylabel("Individual Fibre Transmission") + for fibre_trans in self.output['raw_integration_transmissions']: + ax.plot(self.output['wavelengths'],fibre_trans,color='black',linewidth=0.5) + ax2=ax.twinx() + ax2.axhline(0,label="Individual Fibre",color='black',linewidth=0.5) + ax2.plot(self.output['wavelengths'],np.sum(self.output['raw_integration_transmissions'],axis=0),label="Bundle = {}um".format(self.input['aperture_waveref']),color='red') + plt.axvline(self.input['guide_waveref'],label="Guide = {}um".format(self.input['guide_waveref']),color='black',linestyle='--',linewidth=0.5) + ax2.set_ylabel("Bundle Transmission") + ax2.set_ylim(0,1) + ax.set_xlabel("Wavelength (um)") + plt.legend() + + fig, ax = plt.subplots(figsize=[5,5]) + weights = np.linspace(0, len(self.output['wavelengths'])-1,4) + norm = mpl.colors.Normalize(vmin=min(weights), vmax=max(weights)) + cmap = mpl.cm.ScalarMappable(norm=norm, cmap='seismic') + circle1 = plt.Circle((0, 0), self.output['major_axis']/2, color='black', fill=False, label='~Aperture') + ax.add_patch(circle1) + plt.axvline(0,color='black',linestyle='--',linewidth=0.7,label="PA = {}".format(self.config['relative_plate_PA_angle'])) + plt.scatter(self.output['centre_shift'][0],self.output['centre_shift'][1],label='Guide = {}um'.format(self.input['guide_waveref']),color='black',marker='+') + plt.xlim(-0.4,0.4) + plt.ylim(-0.4,0.4) + shifts=self.output['shifts'] + for i in weights: + xs,ys=[],[] + for o in range(0,len(shifts)): + xs.append(shifts[o][int(i)][0]) + ys.append(shifts[o][int(i)][1]) + plt.plot(xs,ys,marker='x',color=cmap.to_rgba(int(i)),label="{}um".format(round(self.output['wavelengths'][int(i)],4))) + plt.legend() + plt.xlabel("x (arcsec)") + plt.ylabel("y (arcsec)") + + def run(self, HA_start, HA_end, declination, band, + sampling=0.01, guide_waveref=-1, aperture_waveref=-1): + """ + Function to run MOSAIC AD simulation for individual fibre transmission curves. + Goes through multiple steps: + 1) Finds observation airmasses based on provided hour angles and declination + 2) Generates wavelengths to calculate transmission for based on the chosen MOSAIC observing band + 3) Calculates the shifts of these wavelengths at each airmass + 4) Generates aperture (for each fibre) + 5) Generates PSFs of the shifted moffats + 6) Calculates transmissions using apertures and PSFs + + Parameters + ---------- + HA_start,HA_end: float, hours + HA values to start and end the simulated observation at + + declination: float, degrees + Declination of the target in the simulated observation + + band: string + Which MOSAIC observing band to carry out the simulation in. + Must be in the form "X_X_X", such as "LR_VIS_B" with the following options: + {LR: {VIS: {B, G, R, All}, NIR: {IY, J, H, All}}, HR: {VIS: {G, R}, NIR: {IY, H}}} + + sampling: float, um + What interval in um to sample the wavelengths at, default 0.01um + + guide_waveref, aperture_waveref: float, um + What wavelength to guide the telescope on and centre the apertures on halfway through the integration respectively + defaults (when values set to -1) are wavelengths halfway through the chosen band + + Returns + ------- + wavelengths: 1D array of floats + [um] wavelength samples the transmissions have been calculated for + + fibre_transmissions_dic: dictionary, containing 1D arrays of floats + [0..1] Transmissions of each of the fibres for each different + simulated wavelength. Dictionary key is the fibre id, {0,1...}. + The 1D array axis is the wavelength + """ + + self.load_HA(HA_start,HA_end,declination) + self.load_MOSAIC_band(band,sampling) + + if guide_waveref==-1 or aperture_waveref ==-1: + band_mid=round((self.output['wavelengths'][0]+self.output['wavelengths'][-1])/2,3) + guide_waveref,aperture_waveref=band_mid,band_mid + self.input['guide_waveref']=guide_waveref + self.input['aperture_waveref']=aperture_waveref + + self.calculate_shifts(guide_waveref,aperture_waveref) + self.load_aperture(band) + self.load_PSFs() + self.calculate_transmissions() + + wavelengths=self.output['wavelengths'] + fibre_transmissions=self.output['raw_integration_transmissions'] + + fibre_transmissions_dict={} + for count,fibre_transmission in enumerate(fibre_transmissions): + fibre_transmissions_dict[count]=fibre_transmission + + return wavelengths, fibre_transmissions_dict diff --git a/scopesim/effects/mosaic_atmo_disp_fibre_coupling_utils.py b/scopesim/effects/mosaic_atmo_disp_fibre_coupling_utils.py new file mode 100644 index 00000000..f82d43fd --- /dev/null +++ b/scopesim/effects/mosaic_atmo_disp_fibre_coupling_utils.py @@ -0,0 +1,220 @@ +import math +import numpy as np +from astropy.modeling.models import Moffat2D +from astropy.coordinates import Angle +from astropy import units as u +from matplotlib.path import Path +import pandas as pd + + +def calculate_FWHM(wavelength,airmass,config): + median_seeing=config['median_seeing'] + median_seeing_wl=config['median_seeing_wl'] + L0=config['wavefront_outerscale'] + D=config['telescope_diameter'] + + r0=0.1*median_seeing**(-1)*(wavelength/median_seeing_wl)**1.2*airmass**(-0.6) + F_kolb=1/(1+300*(D/L0))-1 + + FWHM_atm=median_seeing*airmass**(0.6)*(wavelength/median_seeing_wl)**(-0.2)*np.sqrt(1+F_kolb*2.183*(r0/L0)**0.356) + FWHM_dl=0.212*wavelength/D + + FWHM_total=np.sqrt(FWHM_atm**2+FWHM_dl**2) + + return FWHM_total + +def make_moffat_PSFs(wavelengths,offsets,airmass,diameter,config,beta=2.5): + scale=config['sim_scale'] + + wavelengths,offsets = np.array(wavelengths),np.array(offsets) + pixel_radius=math.ceil(diameter/2/scale) #radius of aperture in pixels + + FWHMs = calculate_FWHM(wavelengths,airmass,config) + + x = np.arange(0, pixel_radius*2) + y = np.arange(0, pixel_radius*2) + x, y = np.meshgrid(x, y) + + PSFs=np.zeros((len(wavelengths),pixel_radius*2,pixel_radius*2)) + for count in range(0,len(wavelengths)): + alpha=FWHMs[count]/scale/(2*np.sqrt(2**(1/beta)-1)) + moffat_total=(np.pi*alpha**2)/(beta-1) + x_pos=offsets[count][0]/scale + y_pos=offsets[count][1]/scale + + Moffat=Moffat2D(1,x_pos+pixel_radius-0.5,y_pos+pixel_radius-0.5,alpha,beta) + Moffat_data=Moffat(x,y) + PSFs[count]=Moffat_data/moffat_total + + return PSFs + +def parallatic_angle(HA,dec,lat): + HA=Angle(HA*u.hour).rad + q = np.arctan2(np.sin(HA),(np.cos(dec)*np.tan(lat)-np.sin(dec)*np.cos(HA))) + return q + +def diff_shift(wave, airmass, atm_ref_wav,config): + T = config['temperature'] + HR = config['humidity'] + P = config['pressure'] + Lambda0 = atm_ref_wav + wave = wave + + ZD_deg = airmass_to_zenith_dist(airmass) + ZD = np.deg2rad(ZD_deg) + + # saturation pressure Ps (millibars) + PS = -10474.0 + 116.43*T - 0.43284*T**2 + 0.00053840*T**3 + + # water vapour pressure + Pw = HR * PS + # dry air pressure + Pa = P - Pw + + #dry air density + Da = (1 + Pa * (57.90*1.0e-8 - 0.0009325/T + 0.25844/T**2)) * Pa/T + #water vapour density + Dw = (1 + Pw * (1 + 3.7 * 1E-4 * Pw) * (- 2.37321 * 1E-3 + 2.23366/T + - 710.792/T**2 + + 77514.1/T**3)) * Pw/T + + S0 = 1.0/Lambda0 + S = 1.0/wave + + N0_1 = (1.0E-8*((2371.34+683939.7/(130.0-S0**2)+4547.3/(38.9-S0**2))*Da + + (6487.31+58.058*S0**2-0.71150*S0**4+0.08851*S0**6)*Dw)) + N_1 = 1.0E-8*((2371.34+683939.7/(130.0-S**2)+4547.3/(38.9-S**2))*Da + + (6487.31+58.058*S**2-0.71150*S**4+0.08851*S**6)*Dw) + + DR = np.tan(ZD)*(N0_1-N_1) * u.rad + + return (DR.to(u.arcsec)).value + +def airmass_to_zenith_dist(airmass): + return np.rad2deg(np.arccos(1/airmass)) + +def line(A,B): + m=(A[1]-B[1])/(A[0]-B[0]) + c=A[1]-m*A[0] + return m,c + +def make_aperture_old(band,major_axis,scale): + pixel_radius=math.ceil(major_axis/2/scale) #radius of aperture in pixels + + if band[0:6]=="LR_VIS" or band[0:6]=="LR_NIR" or band[0:6]=="HR_NIR": + sampling = major_axis/3/scale + + triangle_side=sampling*np.sqrt(3)/3 + aperture_centre=[pixel_radius-0.5,pixel_radius-0.5] + + centre_0=aperture_centre + centre_1=[centre_0[0],centre_0[1]+sampling] + centre_2=[centre_0[0]+np.sqrt((triangle_side*3/2)**2+(sampling/2)**2)*np.cos(np.pi/6),centre_0[1]+np.sqrt((triangle_side*3/2)**2+(sampling/2)**2)*np.sin(np.pi/6)] + centre_3=[centre_0[0]+np.sqrt((triangle_side*3/2)**2+(sampling/2)**2)*np.cos(np.pi/6),centre_0[1]-np.sqrt((triangle_side*3/2)**2+(sampling/2)**2)*np.sin(np.pi/6)] + centre_4=[centre_0[0],centre_0[1]-sampling] + centre_5=[centre_0[0]-np.sqrt((triangle_side*3/2)**2+(sampling/2)**2)*np.cos(np.pi/6),centre_0[1]-np.sqrt((triangle_side*3/2)**2+(sampling/2)**2)*np.sin(np.pi/6)] + centre_6=[centre_0[0]-np.sqrt((triangle_side*3/2)**2+(sampling/2)**2)*np.cos(np.pi/6),centre_0[1]+np.sqrt((triangle_side*3/2)**2+(sampling/2)**2)*np.sin(np.pi/6)] + + centres=[centre_0,centre_1,centre_2,centre_3,centre_4,centre_5,centre_6] + + apertures=np.zeros((7,pixel_radius*2,pixel_radius*2)) + apertures_table=[] + for count,centre in enumerate(centres): + P1=[centre[0]+triangle_side*np.cos(np.pi*1/3),centre[1]+triangle_side*np.sin(np.pi/3)] + P2=[centre[0]+triangle_side,centre[1]] + P3=[centre[0]+triangle_side*np.cos(np.pi/3),centre[1]-triangle_side*np.sin(np.pi/3)] + P4=[centre[0]-triangle_side*np.cos(np.pi/3),centre[1]-triangle_side*np.sin(np.pi/3)] + P5=[centre[0]-triangle_side,centre[1]] + P6=[centre[0]-triangle_side*np.cos(np.pi/3),centre[1]+triangle_side*np.sin(np.pi/3)] + + polygon=[P1,P2,P3,P4,P5,P6] + height=pixel_radius*2 + width=pixel_radius*2 + poly_path=Path(polygon) + + x, y = np.mgrid[:height, :width] + coors=np.hstack((x.reshape(-1, 1), y.reshape(-1,1))) + + mask = poly_path.contains_points(coors) + mask=mask.reshape(height, width) + + apertures[count]=mask + apertures_table.append([count,(centre[1]-aperture_centre[1])*scale,(centre[0]-aperture_centre[0])*scale]) + + apertures_table=pd.DataFrame(apertures_table,columns=['fibre_id','dx','dy']) + + return apertures,apertures_table + +def make_aperture(band,major_axis,config): + scale=float(config['sim_scale']) + pixel_radius=math.ceil(major_axis/2/scale) #radius of aperture in pixels + + if config['custom_hexagon_radius'] != 0: + radius=config['custom_hexagon_radius'] + else: + radius=config[band[0:6]+"_hexagon_radius"] + sampling = major_axis/(2*radius-1)/scale + triangle_side=sampling*np.sqrt(3)/3 + aperture_centre=[pixel_radius-0.5,pixel_radius-0.5] + + centres_array=(cube_spiral([0,0],radius-1)) + q_vector=np.array([sampling/2,sampling/np.sqrt(3)*3/2]) + r_vector=np.array([sampling,0]) + + centres_xy_coords=[] + for array_val in centres_array: + xy_vals=np.array(array_val[0]*q_vector+array_val[1]*r_vector) + centres_xy_coords.append(xy_vals) + + apertures=np.zeros((len(centres_array),pixel_radius*2,pixel_radius*2)) + apertures_table=[] + for count,vals in enumerate(centres_xy_coords): + centre=[vals[1]+aperture_centre[0],vals[0]+aperture_centre[1]] + P1=[centre[0]+triangle_side*np.cos(np.pi*1/3),centre[1]+triangle_side*np.sin(np.pi/3)] + P2=[centre[0]+triangle_side,centre[1]] + P3=[centre[0]+triangle_side*np.cos(np.pi/3),centre[1]-triangle_side*np.sin(np.pi/3)] + P4=[centre[0]-triangle_side*np.cos(np.pi/3),centre[1]-triangle_side*np.sin(np.pi/3)] + P5=[centre[0]-triangle_side,centre[1]] + P6=[centre[0]-triangle_side*np.cos(np.pi/3),centre[1]+triangle_side*np.sin(np.pi/3)] + polygon=[P1,P2,P3,P4,P5,P6] + height=pixel_radius*2 + width=pixel_radius*2 + poly_path=Path(polygon) + + x, y = np.mgrid[:height, :width] + coors=np.hstack((x.reshape(-1, 1), y.reshape(-1,1))) + + mask = poly_path.contains_points(coors) + mask=mask.reshape(height, width) + + apertures[count]=mask + apertures_table.append([count,(vals[0])*scale,(vals[1])*scale]) + + apertures_table=pd.DataFrame(apertures_table,columns=['fibre_id','dx','dy']) + return apertures,apertures_table + +def cube_direction(direction): + direction_vectors=[[+1,0],[+1,-1],[0,-1],[-1,0],[-1,+1],[0,+1]] + return direction_vectors[direction] +def cube_add(hex,vec): + return [hex[0]+vec[0],hex[1]+vec[1]] +def cube_neighbor(cube,direction): + return cube_add(cube,cube_direction(direction)) +def cube_scale(hex, factor): + return [hex[0]*factor,hex[1]*factor] +def cube_ring(center, radius): + results = [] + hex = cube_add(center, + cube_scale(cube_direction(4), radius)) + for i in range(0,6): + for j in range(0,radius): + results.append(hex) + hex = cube_neighbor(hex, i) + return results +def cube_spiral(center, radius): + results = [center] + for k in range(1,radius+1): + ring_results=cube_ring(center, k) + for i in ring_results: + results.append(i) + return results \ No newline at end of file diff --git a/scopesim/tests/tests_effects/test_MosaicAtmosDispFibreCoupling.py b/scopesim/tests/tests_effects/test_MosaicAtmosDispFibreCoupling.py new file mode 100644 index 00000000..1d3c1668 --- /dev/null +++ b/scopesim/tests/tests_effects/test_MosaicAtmosDispFibreCoupling.py @@ -0,0 +1,32 @@ +import pytest +from pytest import approx + +import numpy as np +from matplotlib import pyplot as plt +from matplotlib.colors import LogNorm +from astropy import units as u + +from scopesim.effects import MosaicAtmosDispFibreCoupling +from scopesim import rc + +PLOTS = False + + +class TestInit: + def test_initialises_with_nothing(self): + fc = MosaicAtmosDispFibreCoupling() + assert isinstance(fc, MosaicAtmosDispFibreCoupling) + + def test_get_output_from_fibre_coupling_class(self): + """ + Not sure if this is a worthy test, but it shows me what I want to see + """ + + fc = MosaicAtmosDispFibreCoupling() + wave, trans = fc.coupling_fractions.run(0, 0, 0, "LR_NIR_H", 0.01, -1, -1) + for tran in trans: + assert np.all(tran > 0) + + +class TestApplyTo: + pass