#!/usr/bin/env python3
"""
This file is part of the PAFFrontendSim.
The PAFFrontendSim is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License v3 as published by the Free Software Foundation.
The PAFFrontendSim is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
See the GNU General Public License for more details.
You should have received a copy of the GNU General Public License along with the PAFFrontendSim. If not, see https://www.gnu.org/licenses/.
"""
#############################################################################
# ToDo
#############################################################################
# - improve readability, add more comments
#
# - implement atmospheric model correctly
#
# - remove object lists and use optics data classes instead
#
#############################################################################
# import libs
#############################################################################
import numpy as np
import matplotlib.pyplot as plt
from copy import copy, deepcopy
import importlib
import inspect
from pathlib import Path
import LogClass
from Simulation import AtmosphereDataClass
from Simulation.UFO import UFODataClass
from Simulation.UFO.UFOGaussCalc import calcConicData, complete
from Simulation.UFO import OpticsDataClass
from Simulation.UFO.UFO_CCode import _UFOVector
from Simulation import SignalChainDataClass
from Simulation.constants import _SkyDistance, _SkySize, _SkyStep, _TELESCOPE_RASTER, _DEFAULT_WL
from Simulation.constants import _Tskatter, _Tiso, _Tgrd, _Tcbg, _Tatm, _Tloss, _Sig, _kB, _ETA, _C, _SigPol
from Simulation.constants import _Bandwidth, _Impedance, _I0iso
from Simulation.constants import _pythonPath
#############################################################################
# constants
#############################################################################
_CLASSNAME = "Telescopes"
#############################################################################
# classes
#############################################################################
[docs]
class Telescope :
"""
This module represents the telescope and defines the surfaces for the optical physics simulation.
It contains methods to generate the primary and secondary reflectors as well as the sky and spillover regions.
Blockage caused by support structures can be taken into account.
Customs specifications and geometries can be imported as external Python modules.
Optionally an atmosphere (through AtmosphereDataClass) and receiver optics (through OpticsDataClass) can be added to the simulation environment.
This module uses the UFO modules to perform the optical physics simulation.
"""
def __init__ ( self, log = None):
"""
Initialize class.
Parameters:
log (LogClass): Logging Class used by the calling application.
"""
if log == None:
self.log = LogClass.LogClass(_CLASSNAME)
else:
self.log = deepcopy(log) # logging class used for all in and output messages
self.log.addLabel(self)
self.Name = None # telescope name
self.mirrorList = [] # list of the telescope mirrors (UFODataClasses) which are simulated one after the other
self.mirrorName = [] # list of mirror ID-strings
self.spillList = [] # list of spill-over areas to compute (UFODataClasses)
self.spillStep = [] # step of telescope simulation to take to calculate spill (-1 == detector)
self.spillTemp = [] # termination temperature of the spill-over region
self.spillNames = [] # list of spill-over region names (should be 'sky' for sky terminated spill-over, 'ground' for ground terminated spill-over, or 'skatter')
self.Diameter = 0. # telescope diameter [mm]
self.FocalLength = 0. # telescope focal length [mm]
self.weather = 10. # perceptible water vaport (average codition) at the site [mm]; used within the analyzing part
self.Elevation = 90. # telescope elevation position [deg] (for the temperature distribution in the spill-over area). Zenith = 90.
self.OpeningAngle = 0. # opening angle [deg] (from telescope edge to focal point)
self.secDiam = None # secondary diameter (None if no secondary is given)
self.RMS = 0. # surface RMS
self.EdgeTaper = -12. # optimum (or to be used Edge taper
self.waist = None # function to determin waist size for optimm illumination , takes wavelength as argument
self.focus = None # function giving the focus position for a particular frequency
self.optics = None # OpticsDataClass describing a fore optics in front of the telescope. Is simulated automatically if present
self.inputZ = ( 0., 0., -1. ) # telescope input coordinate system z-axis (light direction)
self.inputX = ( -1., 0., 0. ) # telescope input coordinate system x-axis (polarization 0 degree position)
self.inputPos = ( 0., 0., 0. ) # telescope input port location (coordinate system base point)
self.atmosphere = None # class for calculating atmospheric transmission and sky temperature (only implemented if not set to None)
self.usedFolder = "" # folder for debug-file writing, no file writing if empty
self.Tsky = _Tcbg + _Tatm
[docs]
def printOut ( self, line, message_type, prefix = True, lineFeed = True ):
"""
Wrapper for the printOut method of the actual log class.
"""
self.log.printOut(line, message_type, prefix, lineFeed)
[docs]
def importAtmosphere( self, folder, name, season, pwv ):
"""
Attempts to import new atmospheric data for using AtmosphereDataClass.
If the file doesn't exist then nothing is done.
Parameters:
folder (string): Path of the folder the file is expected in
name (string): Name of the telescope
season (string): Season of the data (in 3 letter format)
pwv (string): Perceptible water vapor in percent
"""
targetFilePath = Path(folder, name+"_"+season+"_"+pwv+".out") # Constructing string of full path
if targetFilePath.exists(): # If the file exists
self.atmosphere = AtmosphereDataClass.Atmosphere( log = self.log )
self.atmosphere.readFile(targetFilePath) # Read it
self.Tsky = 0.
else: # Else file doesn't exist
self.printOut("File not found: %s" % targetFilePath, self.log.ERROR_TYPE)
[docs]
def importTelescopeDict(self, moduleString, name):
"""
Imports a Python module containing the definition of a telescope.
Modules must define a telescope either through a dictionary of specifications or a function.
Parameters:
moduleString (string): Import string of the module.
name (string): Name of the telescope.
Returns:
telescopeDataDict (dict): Contains either the specifications of the telescope or a description + function call.
"""
imported_module = importlib.import_module(moduleString) # Import module
for obj_name, obj in inspect.getmembers( imported_module ) : # loop over all members of the module
if obj_name == name : # check if member name is identical to file-name ==> dictionary of specifications
moduleType = moduleString.split(".")[-2]
if moduleType == "Specifications":
telescopeDataDict = obj # Store the dictionary of specifications
elif moduleType == "Functions":
telescopeDataDict = {
"Description" : imported_module._description_, # Store the description within the module
"FunctionCall" : deepcopy( obj ), # Store the function call
} # Store description and function call in a dictionary
else:
self.printOut("Telescopes must be defined either through 'Functions' or 'Specifications', but got '%s'" % moduleString.split(".")[-2], self.log.ERROR_TYPE)
telescopeDataDict = None
return telescopeDataDict
[docs]
def attemptTelescopeDataImport(self, telescopePath, moduleType, name):
"""
Attempts to retrieve a module of telescope data.
It checks if the file exists and uses self.importTelescopeDict() to import the data in said file.
Parameters:
telescopePath (Path): The path preceding the folder in which the module is expected to be located.
moduleType (String): Name of the module's folder. Must be either "Specifications" or "Functions".
name (String): Name of the module, excluding the ".py" file-extension.
Returns:
telescopeDataDict (dict): Contains the data found in the module, otherwise None.
"""
filePath = Path(telescopePath, Path(moduleType, name+".py")) # reassemble the full path of the module
if filePath.exists(): # if the file exists
moduleString = telescopePath.parts[-2] + "." + telescopePath.parts[-1] + "." + moduleType + "." + name # assemble the import string of the module
telescopeDataDict = self.importTelescopeDict(moduleString, name) # Store the specification dictionary with that name
self.printOut("Found and imported %s" % Path(*filePath.parts[-5:]), self.log.FLOW_TYPE)
else:
telescopeDataDict = None
self.printOut("Couldn't find %s" % Path(*filePath.parts[-5:]) , self.log.FLOW_TYPE)
return telescopeDataDict
[docs]
def createWithFunction(self, wavelength, CalcSpill, RMS, spillPrim):
"""
Creates a telescope using an external function.
If a corresponding specification dictionary exists then it is passed to the function.
The function might require it for calculations and it might return the dictionary with modified specifications.
Parameters:
wavelength (float): spacial rastering of the telescope along this dimension
CalcSpill (boolean): Add spill over calculation to the telescope model (only useful if a noise model is established)
RMS (float): Main dish RMS (if negative telescope nominal RMS is used)
spillPrim (boolean): indicating if primary spill-over area is used or not
"""
if self.specsDict != None: # If a corresponding dictionary of specifications exists
self.specsDict["Description"] = self.functionDict["Description"] # Overwrite its description
if ( self.specsDict.get( 'CreateWithFunc', True ) ) : # check if telescope should be created according to specs before function call
# otherwise the function is expected to generate the entire telescope itself
self.createFromSpecs(self.specsDict, wavelength, CalcSpill, RMS, spillPrim) # Create the telescope based on specs
self.functionDict["FunctionCall"]( self, wavelength, CalcSpill, RMS, spillPrim, specs = self.specsDict ) # Run function of external module and store returned specs
[docs]
def createFromSpecs(self, telDict, wavelength, CalcSpill, RMS, spillPrim ):
"""
Create the telescope assembly.
Parameters:
wavelength (float): spacial rastering of the telescope along this dimension
CalcSpill (boolean): Add spill over calculation to the telescope model (only useful if a noise model is established)
RMS (float): Main dish RMS (if negative telescope nominal RMS is used)
spillPrim (boolean): indicating if primary spill-over area is used or not
"""
self.printOut("description : %s" % ( telDict["Description"] ), self.log.DATA_TYPE, False)
if 'TelTaper' in telDict['kargs']:
self.printOut("Edge-taper : %.1f dB" % ( telDict['kargs']['TelTaper'] ), self.log.DATA_TYPE, False )
# Perform the corresponding telescope creation function
if telDict["Module"] == "GregorianOrCass" : # Gregorian or Cassegrain with conic surfaces (secondary focus)
self.GregorianOrCass( wavelength, CalcSpill, RMS, spillPrim, **telDict["kargs"] )
elif telDict["Module"] == "Parabolic" : # Parabolic main dish (primary focus)
self.Parabolic ( wavelength, CalcSpill, RMS, spillPrim, **telDict["kargs"] )
else:
self.printOut("Tried to create a telescope of unknown type > %s <" % telDict["Module"] , self.log.ERROR_TYPE)
[docs]
def setTelescope( self, name, wavelength = _DEFAULT_WL, folder = "", CalcSpill = True, Raster = 1, RMS = 0, spillPrim = True, Elevation = 90. ):
"""
Set and create a telescope
Parameters:
name (string): The name of the telescope we want to use
wavelength (float): spacial rastering of the telescope along this dimension
folder (string): If given, folder where the telescope data files are stored for debugging
CalcSpill (boolean): Add spill over calculation to the telescope model (only useful if a noise model is established)
Raster (float): Multiplicator for internal telescope rastering (sometimes required for higher frequencies)
RMS (float): Main dish RMS (if negative telescope nominal RMS is used)
spillPrim (boolean): indicating if primary spill-over area is used or not
Elevation (float): Elevation in degrees at which the telescope is created. 90deg corresponds to zenith (default)
"""
file_path = Path(__file__).resolve() # get the path of this file
full_path_tuple = file_path.parts # split the path into its components
root_index = full_path_tuple.index("paffrontendsim-d31") # determine at what position in the path the root folder is
root_path_tuple = full_path_tuple[:root_index+1] # store the path of the root folder
telescopePath = Path(*root_path_tuple, Path("src", "Libs", "Simulation", "Telescope"))
self.specsDict = self.attemptTelescopeDataImport(telescopePath, "Specifications", name )
self.functionDict = self.attemptTelescopeDataImport(telescopePath, "Functions", name )
if ( self.specsDict is None ) and ( self.functionDict is None ):
self.printOut("Couldn't import any data for %s" % name, self.log.ERROR_TYPE)
self.usedFolder = folder
# delete previous data
self.mirrorList = [] # list of the telescope mirrors (UFODataClasses) which are simulated one after the other
self.mirrorName = [] # list of mirror ID-strings
self.spillList = [] # list of spill-over areas to compute (UFODataClasses)
self.spillStep = [] # step of telescope simulation to take to calculate spill (-1 == detector)
self.spillTemp = [] # termination temperature of the spill-over region
self.spillNames = [] # list of spill-over region names (should be 'sky' for sky terminated spill-over, 'ground' for ground terminated spill-over, or 'skatter')
self.resultList = [] # simulation results
self.secDiam = None # secondary mirror existance flag (diameter if not None)
self.Raster = _TELESCOPE_RASTER * Raster
self.Name = name # Store the name of the telescope we want to generate
self.Elevation = Elevation # set telescope elevation
self.printOut("Setting up telescope %s:" % self.Name, self.log.FLOW_TYPE)
self.printOut("wavelength : %.2f mm" % ( wavelength ), self.log.DATA_TYPE, False)
self.printOut("Surface-RMS : %.3f mm" % ( RMS ), self.log.DATA_TYPE, False)
self.printOut("using spill-over in general: %s" % ( CalcSpill ), self.log.DATA_TYPE, False)
self.printOut("using primary spill-over : %s" % ( spillPrim ), self.log.DATA_TYPE, False)
self.printOut("folder : %s" % ( folder ), self.log.DATA_TYPE, False)
if self.functionDict != None: # If there is an external function
self.printOut("Creating telescope > %s < using function" % ( name ), self.log.FLOW_TYPE, False)
self.createWithFunction(wavelength, CalcSpill, RMS, spillPrim) # Create the telescope using the external function
elif self.specsDict != None: # Else if a specifiation dictionary exists
self.printOut("Creating telescope > %s < based on specs" % ( name ), self.log.FLOW_TYPE, False)
self.createFromSpecs(self.specsDict, wavelength, CalcSpill, RMS, spillPrim) # Create telescope based on specs
else: # Else there are no specs or modules with that name, error
self.printOut("No modules for telescope > %s < found!" % ( name ), self.log.ERROR_TYPE, False)
#####################################
# plotting
[docs]
def plotTelescope ( self, figname = "", show = False ) :
"""
Plot the telescope projections.
Parameters:
figname (string): file-name of the figure (if given figure is only shown on screen if show == True
show (boolean): flag to show plot on screen even when being written to disk (True)
"""
self.printOut("Plotting telescope projection", self.log.FLOW_TYPE)
# setup plot
fig, ax = plt.subplots(2)
ax[0].set_xlabel("distance [mm]")
ax[0].set_ylabel("distance [mm]")
ax[0].set_title ("Projected active telescope area (Top: x-y plane, Bottom, x-z plane)")
ax[0].set_aspect('equal')
ax[1].set_xlabel("distance [mm]")
ax[1].set_ylabel("distance [mm]")
ax[1].set_aspect('equal')
colors = ('r.', 'b.', 'g.', 'k.') # colors of the individual mirrors of the telescope
# plot all mirrors
for obj_idx, obj in enumerate(self.mirrorList[:-1]) : # loop all telescope objects except the sky
x = []
y = []
z = []
for idx in range(0, obj.points) :
x.append(obj.Data[idx].v.x)
y.append(obj.Data[idx].v.y)
z.append(obj.Data[idx].v.z)
x = np.asarray(x)
y = np.asarray(y)
z = np.asarray(z)
ax[0].plot(x, y, colors[obj_idx % len(colors)], label=self.mirrorName[obj_idx])
ax[1].plot(x, z, colors[obj_idx % len(colors)])#, label=self.mirrorName[obj_idx])
dist = self.mirrorList[0].Header.v.z
x = [ -np.sin(self.OpeningAngle/180*np.pi)*dist, 0., np.sin(self.OpeningAngle/180*np.pi)*dist ]
y = [ 0., 0., 0.]
z = [ np.cos(self.OpeningAngle/180*np.pi)*dist, 0., np.cos(self.OpeningAngle/180*np.pi)* dist]
ax[0].plot(x, y, 'y-', label="Opening angle")
ax[0].plot(y, x, 'y-')
ax[1].plot(x, z, 'y-')
ax[1].plot(y, z, 'y-')
fig.legend(loc='upper left', bbox_to_anchor=(0.65, 0.85))
if figname != "" :
fig.savefig(figname)
if show :
plt.show()
else :
plt.show()
plt.close()
#####################################
# using the telescope and data analysis
[docs]
def addOptics ( self, optics = None, FocalPoint = [0, 0, 0] ) :
"""
Adds a fore-optics to the telescope.
Parameters:
optics (OpticsData): optics to be added. The output port of the optics is connected to the telescope, simulation starts at the input port. If no optics is given, an existing optics is removed.
FocalPoint (3d point, float): point in the focal plane at which the receiver beam is entering (usually the center, but for of-axis optics an offset can be given).
"""
if not isinstance(optics, OpticsDataClass.OpticsData):
self.printOut("Deleting optics from telescope", self.log.FLOW_TYPE)
self.optics = None
return
if isinstance(FocalPoint, (int, float)) :
FocalPoint = [FocalPoint]
while len(FocalPoint) < 3 :
FocalPoint.append(0)
FocalPoint = np.asarray(FocalPoint)
self.printOut("Adding optics to telescope", self.log.FLOW_TYPE)
# set optics and move it to telescope input port
self.optics = optics
# get distance to first telescope mirror
distFirstMirror = self.mirrorList[0].Header.v.z
aX = np.arctan2( -FocalPoint[1], -distFirstMirror )
aY = np.arctan2( FocalPoint[0], -distFirstMirror )
rx = np.asarray([
[1, 0, 0],
[0, np.cos(aX), -np.sin(aX)],
[0, np.sin(aX), np.cos(aX)]
])
ry = np.asarray([
[ np.cos(aY), 0, np.sin(aY)],
[ 0, 1, 0],
[-np.sin(aY), 0, np.cos(aY)]
])
outputZ = np.dot( ry, np.dot( rx, np.asarray( ( 0., 0., -1. ) ) ) )
outputX = np.dot( ry, np.dot( rx, np.asarray( ( -1., 0., 0. ) ) ) )
self.optics.setOutputPort( outputZ = outputZ, outputX = outputX, outputPos = FocalPoint )
#set input port to optics input
self.inputZ = self.optics.InputZ
self.inputX = self.optics.InputZ
self.inputPos = self.optics.InputPosition
[docs]
def simulate ( self, farfield, scaleWaveLength = -1 ) :
"""
Simulate the far-field pattern 'farfield' and return sky and spill-over region data.
Parameters:
farfield (FarFieldDataClass): With the farfield of the antenna
scaleWaveLength (int): Wavelength to scale the sky-size with (if negative, farfield wavelength is taken)
Returns:
sky (UFODataClass) : with type 'CURRENT' of the sky beam-pattern
spill list of (UFODataClass) : with type 'CURRENT' of the spill-over area without temperature weighting
spillT list of (UFODataClass) : with type 'CURRENT' of the spill-over area with temperature weighting
"""
self.printOut("Simulating far-field pattern on telescope %s" % self.Name, self.log.FLOW_TYPE)
# check if mirrors are available
if len(self.mirrorList) == 0 :
self.printOut("no Telescope set, exiting", self.log.ERROR_TYPE, False)
return None, None, None
# check if optics is present and if so simulate it
if isinstance(self.optics, OpticsDataClass.OpticsData ) :
self.printOut("using fore-optics", self.log.FLOW_TYPE, False)
opticsOut, opticsSpill, opticsSpillT = self.optics.simulateOptics ( Source = farfield )
ufo = True
farfield = opticsOut
ufo = False
if isinstance(farfield, UFODataClass.UFOData):
ufo = True
if farfield.Wavelength == None :
if scaleWaveLength < 0:
self.printOut("using UFO-Data, but no wavelength is given! Exiting", self.log.ERROR_TYPE, False)
return None, None, None
farfield.Wavelength = scaleWavelength
self.printOut("using UFO-Data.", self.log.DEBUG_TYPE, False)
# deepcopy mirror-list to be able to re-use
skyOrg = deepcopy(self.mirrorList[-1])
sWavelength = scaleWaveLength # TODO too many "wavelength" variables! farfield.Wavelength, sWavelength, scaleWaveLength
if sWavelength < 0 :
sWavelength = farfield.Wavelength
# scale sky (last mirror in mirror-list):
v = self.mirrorList[-1].getCenter()
for idx in range(0, self.mirrorList[-1].points) :
self.mirrorList[-1].Data[idx].v.x = (self.mirrorList[-1].Data[idx].v.x - v[0])*sWavelength/_DEFAULT_WL + v[0]
self.mirrorList[-1].Data[idx].v.y = (self.mirrorList[-1].Data[idx].v.y - v[1])*sWavelength/_DEFAULT_WL + v[1]
self.mirrorList[-1].Data[idx].v.z = (self.mirrorList[-1].Data[idx].v.z - v[2])*sWavelength/_DEFAULT_WL + v[2]
self.mirrorList[-1].Data[idx].df *= sWavelength**2/_DEFAULT_WL**2
######################
# simulate signal path
farfield.setLog( self.log ) # updata farfield log class
resultList = [] # list with all resulting current distributions
self.printOut("simulation of active mirrors :", self.log.FLOW_TYPE, False)
if not ufo : # check for input data type
objCurrent = farfield.calculateCurrent(self.mirrorList[0]) # simulate far-field pattern to first mirror
else :
objCurrent = farfield.simulate( self.mirrorList[0], -1) # simulate ufo-current distribution to first mirror (using internal wavelength)
if self.usedFolder != "" :
objCurrent.writeFile(Path(self.usedFolder, "telCurrent.tclass"))
resultList.append(objCurrent) # add current distribution to result-list
for obj_idx, obj in enumerate(self.mirrorList[1:]) : # loop all mirrors
self.printOut("-- simulating object '%s'"% self.mirrorName[obj_idx+1], self.log.DATA_TYPE, False)
objCurrent = objCurrent.simulate(obj, farfield.Wavelength) # simulate next mirror
resultList.append(objCurrent)
if self.usedFolder != "" :
objCurrent.writeFile(Path(self.usedFolder, "TelescopeMirror_%d.curr" % (obj_idx+1)) )
# rotate sky back in x-y plane
resultList[-1].resetPosition()
# check for atmospheric model
spillTemp = deepcopy( self.spillTemp )
if isinstance( self.atmosphere, AtmosphereDataClass.Atmosphere ) :
atmTrans, tsky = self.atmosphere.SkyTransmission( _C / farfield.Wavelength * 1e-6 ) # C: m/s, wavelength in mm result: GHz
# and adapt sky temperature
for sT in spillTemp:
if sT == self.Tsky :
sT = tsky
lossAtm, lossAtmT = resultList[-1].applyLosses(lossFactor = 1. - atmTrans, lossTemperature = tsky)
else :
lossAtm = UFODataClass.UFOData( log = self.log )
lossAtmT = UFODataClass.UFOData( log = self.log )
lossAtm.Name = 'atmosphere'
lossAtmT.Name = 'atmosphere'
######################
# simulate spill-over
spillList = deepcopy(self.spillList)
# simulate all spill-over areas
self.printOut("simulation spill-over regions :", self.log.FLOW_TYPE, False)
resultSpill = [] # list with all resulting spill-over current distributions
for s_idx, sp in enumerate(spillList) : # loop all spill-over regions
if self.spillStep[s_idx] < 0 : # check if using far-field
self.printOut("-- simulating from detector", self.log.DATA_TYPE, False)
if not ufo: # check for input data type
spillCurrent = farfield.calculateCurrent(sp) # simulate from far-field
else :
spillCurrent = farfield.simulate(sp, -1) # simulate ufo-current distribution to first mirror (using internal wavelength)
resultSpill.append(spillCurrent)
else : # or from current-distribution
self.printOut("-- simulating from element '%s'" % self.mirrorName[self.spillStep[s_idx]], self.log.DATA_TYPE, False)
spillCurrent = deepcopy(resultList[self.spillStep[s_idx]]) # simulate from telescope surface
spillCurrent = spillCurrent.simulate(sp, farfield.Wavelength)
resultSpill.append(spillCurrent)
if self.usedFolder != "" :
spillCurrent.writeFile(Path(self.usedFolder, "Spill-Over_%d.curr" % ( s_idx )))
# combine all spill-over areas
self.printOut("combining spill-areas with temperatures", self.log.DEBUG_TYPE, False)
if len(resultSpill) > 0 :
spill = deepcopy(resultSpill[0])
else:
spill = UFODataClass.UFOData( log = self.log )
spill.Data = []
spill.Name = 'telescope'
spillT = deepcopy(spill)
for sp_idx, sp in enumerate(resultSpill) : # loop all spill-areas
if self.spillTemp[sp_idx] > 0: # ignore intermediate results
spT = deepcopy(sp)
spT._modField ( np.sqrt( self.spillTemp[sp_idx] / _Tiso ), onlyE = False )
for idx in range(0, sp.points) : # loop data
spill.Data.append(sp.Data[idx])
spillT.Data.append(spT.Data[idx])
# calculate missing energy (e.g from scattering to sky outside the FoV of the sky-plane)
# ToDo
# set ufo-data classes
spill.points = len(spill.Data)
spill.Data = (UFODataClass.UFODataPoint * spill.points)(*spill.Data)
spillT.points = len(spillT.Data)
spillT.Data = (UFODataClass.UFODataPoint * spillT.points)(*spillT.Data)
self.resultList = resultList
self.mirrorList[-1] = skyOrg
#return resultList[-1], spill, spillT
return resultList[-1], [spill, lossAtm], [spillT, lossAtmT]
[docs]
def simulateAndAnalyze ( self, farfield, scaleWavelength = -1, SignalChain = None, center = [0., 0., 0.], useMax = False, lossTemp = -1, lossField = None, lossFieldT = None, inputPower = None) :
"""
Outdated, check if it is in still in use and delete the function
Simulate the far-field pattern 'farfield' and return sky and spill-over region data.
Parameters:
farfield (FarFieldDataClass): With the farfield of the antenna
scaleWavelength (float): Wavelength to scale the sky-size with (if negative, farfield wavelength is taken)
SignalChain (SignalChainDataClass): To calculate rx noise and so on
center (list): 3D vector defining the beam-center for the analysis
useMax (boolean): True indicates to use the maximum intensity as beam-center
lossTemp (float): Temperature [K] at which the losses are realized (negative value uses the default from PAF-simulations)
lossField (ufoDataClass): with field data representing the losses
lossFieldT (ufoDataClass): ufoDataClass with field data representing the temperature correted losses. lossField and lossFieldT must be given together to have an effect
inputPower (float or None): positive value gives the inputPower into the telescope, None: input power is calculated from input surface, negative: use full power as from excitation current
Returns:
tuple (tuple): sky, spill, spillT, DataPack -- UFOData of the simulation, and a list of analysis data
"""
def autoCorrelate ( ufo, center = [0., 0., 0.], useMax = False ):
"""
evaluate the complex auto-correlation function of the UFO data class
"""
center = np.asarray(center)
if useMax:
maxV = 0.
maxIdx = -1
for idx in range(0, ufo.points) :
if ufo.Data[idx].I > maxV :
maxIdx = deepcopy(idx)
maxV = ufo.Data[idx].I
if maxIdx < 0 :
center = [0., 0., 0.]
else:
center = np.asarray([ufo.Data[maxIdx].v.x, ufo.Data[maxIdx].v.y, ufo.Data[maxIdx].v.z])
else :
center = np.asarray([ufo.Header.v.x + center[0], ufo.Header.v.y + center[1], ufo.Header.v.z] + center[2])
rmin = 1e100
Ecenter = None
A = 0. + 0.*1j
for idx in range( 0, ufo.points ) :
E = np.asarray( [ ufo.Data[idx].Er.x+1j*ufo.Data[idx].Ei.x, ufo.Data[idx].Er.y+1j*ufo.Data[idx].Ei.y, ufo.Data[idx].Er.z+1j*ufo.Data[idx].Ei.z ] )
A += np.dot( E, E.conj() ) * ufo.Data[idx].df
r_ = np.asarray( [ ufo.Data[idx].v.x - center[0], ufo.Data[idx].v.y - center[1], ufo.Data[idx].v.z - center[2] ] )
r = np.dot(r_,r_)
if r < rmin :
Ecenter = E
rmin = r
A *= (0.5/_ETA)
return A, Ecenter
# Main part
TlossRx = _Tloss
if lossTemp >= 0. : TlossRx = lossTemp
# simulate telescope
sky, spill, spillT = self.simulate( farfield, scaleWavelength )
sky.resetPosition()
# extract simulation data
if farfield != None :
Wavelength = farfield.Wavelength
else:
farfield = self.optics.Elements[0].ufo
Wavelength = farfield.Wavelength
if Wavelength == None:
Wavelength = scaleWavelength
# get sky transmission and sky-temperature
if isinstance( self.atmosphere, AtmosphereDataClass.Atmosphere ):
skyTrans, skyTemp= self.atmosphere.SkyTransmission( frequency = _C / Wavelength * 10**(-9))#, elevation = self.Elevation, pwv = self.weather)
else :
skyTrans = 1.
skyTemp = self.Tsky
Wavelength = Wavelength / 1000. # wavelength in m not in mm
Aphys = np.pi*self.Diameter**2./4. / 1000000. # effective area im m^2 not mm^2
if isinstance(farfield, UFODataClass.UFOData):
I0 = farfield.I0 / 1000.
if I0 == None :
I0 = np.sqrt(2.*farfield.getPower()[0]/_Impedance) / 1000.# extract excitation current
Impedance = _Impedance
DetectorGain = 1
B = _Bandwidth
if type(inputPower) == type(None):
inputPower = farfield.getPower()[0] / 1000000.
else:
I0 = farfield.I0
Impedance = farfield.Impedance
DetectorGain = farfield.gain
B = farfield.Bandwidth
if type(inputPower) == type(None):
inputPower = farfield.getPower()
refPower = 0.5 * I0**2 * Impedance * 1.e6
if inputPower < 0. :
inputPower = refPower
I0Input = np.sqrt( 2. * inputPower / _Impedance ) # virtual excitation current for an idealized receiver
self.printOut("Analysis:", self.log.DATA_TYPE)
self.printOut("Input-Power : %.2f" % inputPower, self.log.DATA_TYPE, False)
self.printOut("Reference-Power : %.2f" % refPower , self.log.DATA_TYPE, False)
# set up correclation matrices
if isinstance(lossField, UFODataClass.UFOData) and isinstance(lossFieldT, UFODataClass.UFOData):
lossField, E = autoCorrelate(lossField, center, useMax)
lossFieldT, E = autoCorrelate(lossFieldT, center, useMax)
else :
lossField = 0. + 0*1j
lossFieldT = 0. + 0*1j
if sky == None :
return
Asky , Ec = autoCorrelate(sky, center, useMax) # used only to extract the source E-field vector Ec
Aspill, E = autoCorrelate(spill, center, useMax)
AspillT, E = autoCorrelate(spillT, center, useMax)
# correct for scatter loss, sky simulation only covers the main beam with a roughly 10 omega contour. Full sky power must be difference between input power and spill-power
#if inputPower > 0 :
#Asky = inputPower - Aspill
A = ( Asky + Aspill ) / 1000000. # factor 1/10^6: r^2 in in mm not in m
Asky = Asky / 1000000.
Aspill = Aspill / 1000000.
AspillT = AspillT / 1000000.
# calculate square of the radiation efficiency via the input power
radEff2 = A / inputPower
Aloss = ( 1. - radEff2 )
# rx-losses
AlossRec = lossField / 1000000.
AlossRecT = lossFieldT / 1000000.
# Calculate open circuit thermal noise resistances
# Isotropic noise
Riso_oc = 8.* _kB * B * Impedance * A * _Tiso / inputPower # isotropic noise resistance of the signals reaching the telescope
# external noise
RskyT_oc = 8.* _kB * B * Impedance * Asky * skyTemp / inputPower # sky noise resistance
# spill-over
RspillT_oc = 8.* _kB * B * Impedance * AspillT * _Tiso / inputPower # ground spill-over noise resistance (including temperature variations)
Rspill_oc = 8.* _kB * B * Impedance * Aspill * _Tiso / inputPower # spill-over noise resistance ( for efficiencies )
# loss noise (all missing parts of energy compared to input-power)
RlossT_oc = 8.* _kB * B * Impedance * Aloss * skyTemp # this is scattered power, terminating it on sky
Rloss_oc = 8.* _kB * B * Impedance * Aloss * _Tiso # this is scattered power
# receiver noise
RrecLossT_oc = 8.* _kB * B * Impedance * AlossRecT * _Tiso / refPower # receiver loss noise resistance (including temperature variations)
RrecLoss_oc = 8.* _kB * B * Impedance * AlossRec * _Tiso / refPower # receiver loss noise resistance ( for efficiencies )
# thermal noise (contributions from isotropic temperature and scatter-loss mechanisms) Here all noise in respect to reference power is required
Rt_oc = ( Riso_oc + Rloss_oc ) * inputPower / refPower + RrecLoss_oc
# Rx input impedance using resonant approximation (X = 0) and taking all thermal noise contributions in front of the signal chain input into account
ZA = Rt_oc / ( 8. * _kB * _Tiso * B )
ZA_H = np.transpose(ZA.conj()) # Hermetian of ZA
# Receiver noise
# get signal chain data for the actual wavelength
if not isinstance(SignalChain, SignalChainDataClass.SignalChain) :
SignalChain = SignalChainDataClass.SignalChain( log = self.log )
SignalChain.makeDataset ("Universal", [0.5, 10, 18., 26., 34., 50., 67., 116., 230., 345, 460., 690 ], Tmin = [3., 5., 10., 12., 15., 20., 24., 30., 40., 60., 100., 130.])
f = _C / Wavelength * 10**(-9.) # frequency in GHz
tmin, rn, zopt, s11, s12, s21, s22 = SignalChain.getData(f)
VR2 = 4.* _kB * SignalChain.T0 * rn
Yc = tmin /(2.* SignalChain.T0)/ rn
# analog signal chain input impedance matrix
ZS = SignalChain.Zs
# using antenna array self-impedance for matching
Yopt = 1/ZA
Zopt = ZA
# calculate open circuit receiver noise matrix
Yc = Yc - Yopt;
Yc_H = np.transpose(Yc.conj())
IR2 = np.dot(np.absolute(Yopt)**2., VR2);
Rrec_oc = VR2 + np.dot(np.dot(ZA, Yc), VR2) # receiver noise covariance matrix
Rrec_oc += np.dot(np.dot(VR2, Yc_H), ZA_H) # continuing
Rrec_oc += np.dot(np.dot(ZA, IR2), ZA_H) # continuing
Rrec_oc = 2.* B * Rrec_oc # continuing
# calculate signal response matrix including RFI
# Signal response
k = 2.*np.pi/ Wavelength # wavenumber (rad/m)
E_sig = np.sqrt(_ETA * _Sig * 1e-26 * B ) # electric field intensity in one polarization (V/m)
c1 = 4.* np.pi *1j/(k * _ETA * I0Input / 1000.) # reciprocity theorem scale factor
Vnorm = c1*E_sig / DetectorGain / 1e6 # normalization + take out detector gain
exFactor = np.cos(np.pi*_SigPol/180.)
eyFactor = np.sin(np.pi*_SigPol/180.)
ds = np.dot(Ec, np.array([exFactor, eyFactor, 0]))*_SkyDistance
Vsig_oc = Vnorm * np.conj(ds) # open circuit voltage; simulation adds fields, therefore phase is shifted by 180deg to receiving
# signal correlation
Rsig_oc = np.conj(Vsig_oc).T * Vsig_oc * skyTrans # signal
# Transform all open circuit voltages to loaded array output voltages
# transfer matrix (antenna voltage output -> signal chain input)
Q = ZS/(ZS + ZA) # transfer from antenna input to signal chain input
QH = np.transpose(Q.conj()) # Hermetian of Q
# Transfer all correlation matrices from open circuit to LNA input
Rsig = np.dot(np.dot(Q, Rsig_oc), QH)
Vsig = np.dot(Q, Vsig_oc.conj().T)
RskyT = np.dot(np.dot(Q, RskyT_oc), QH)
Rspill = np.dot(np.dot(Q, Rspill_oc), QH)
RspillT = np.dot(np.dot(Q, RspillT_oc), QH)
Rloss = np.dot(np.dot(Q, Rloss_oc), QH)
RlossT = np.dot(np.dot(Q, RlossT_oc), QH)
Rrec = np.dot(np.dot(Q, Rrec_oc), QH)
RrecLoss = np.dot(np.dot(Q, RrecLoss_oc), QH)
RrecLossT = np.dot(np.dot(Q, RrecLossT_oc), QH)
Rt = np.dot(np.dot(Q, Rt_oc), QH)
Riso = np.dot(np.dot(Q, Riso_oc), QH)
Rnoise = RskyT + RspillT + RlossT + Rrec + RrecLossT
self.printOut("Noise resistances:", self.log.DATA_TYPE)
self.printOut("Rsig : %.2f" % Rsig, self.log.DATA_TYPE, False)
self.printOut("RskyT : %.2f" % RskyT, self.log.DATA_TYPE, False)
self.printOut("RspillT : %.2f" % RspillT, self.log.DATA_TYPE, False)
self.printOut("RlossT : %.2f" % RlossT, self.log.DATA_TYPE, False)
self.printOut("Rrec : %.2f" % Rrec, self.log.DATA_TYPE, False)
self.printOut("RrecLossT : %.2f" % RrecLossT, self.log.DATA_TYPE, False)
self.printOut("Rnoise : %.2f" % Rnoise, self.log.DATA_TYPE, False)
self.printOut("Rt : %.2f" % Rt, self.log.DATA_TYPE, False)
self.printOut("Riso : %.2f" % Riso, self.log.DATA_TYPE, False)
# calculate efficiencies
eta_sp = 1. - np.absolute( np.real( Rspill / Riso ) ) # spill-over efficiency (Riso is the full power reaching the telescope)
eta_rad = np.absolute(np.real( Riso / Rt))
eta_ae = np.absolute( np.real( _Tiso * Rsig / Riso ) )
eta_ae = eta_ae * _kB / ( 0.5 * Aphys * _Sig * 1e-26 )
# calculate temperatures from noise matrices
Tsys = np.absolute(np.real(_Tiso * Rnoise / Rt))
Trec = np.absolute(np.real(_Tiso * Rrec / Rt))
TrecLoss = np.absolute(np.real(_Tiso * RrecLossT / Rt))
Tsp = np.absolute(np.real(_Tiso * RspillT / Rt))
Tsky = np.absolute(np.real(_Tiso * RskyT / Rt))
Tloss = np.absolute(np.real(_Tiso * RlossT / Rt))
Tsig = np.absolute(np.real(_Tiso * Rsig / Rt))
Tsys_oe = Tsys/eta_ae
DataPack = [
eta_ae, eta_sp, eta_rad, Tsys_oe, Tsys, Trec + TrecLoss, Trec, Tsp, Tloss, Tsky, Tsig, (A + Aloss + AlossRec)/refPower, Asky/refPower, Aspill/refPower, Aloss/refPower, AlossRec/refPower
]
self.printOut("eta ae : %.2f" % eta_ae, self.log.DATA_TYPE, False)
self.printOut("eta spill : %.2f" % eta_sp, self.log.DATA_TYPE, False)
self.printOut("eta rad : %.2f" % eta_rad, self.log.DATA_TYPE, False)
self.printOut("Trec-loss : %.2f K" % TrecLoss, self.log.DATA_TYPE, False)
self.printOut("Trec : %.2f K" % Trec, self.log.DATA_TYPE, False)
self.printOut("Tsys : %.2f K" % Tsys, self.log.DATA_TYPE, False)
self.printOut("Tspill : %.2f K" % Tsp, self.log.DATA_TYPE, False)
self.printOut("Tloss : %.2f K" % Tloss, self.log.DATA_TYPE, False)
self.printOut("Tsky : %.2f K" % Tsky, self.log.DATA_TYPE, False)
self.printOut("Tsig : %.2f K" % Tsig, self.log.DATA_TYPE, False)
self.printOut("Tsys/eta : %.2f K" % Tsys_oe, self.log.DATA_TYPE, False)
return sky, spill, spillT, DataPack
[docs]
def makeSpill ( self, ufo, center = [0., 0., 0.] ) :
"""
Points given in the ufo-object are arranged on a sphere with the center center.
The system is rotated until the refPoint and the sphere center are aligned on the z-axis.
All z coordinates are then replaced by values bringing them on the shpere.
Sphere radius is the length of the difference vector between refPoint and center.
Parameters:
ufo (ufo-object): to modify
center (list): vector pointing to the spill-over center
"""
newR = _SkyDistance
if ufo.points > 0 :
center = np.asarray(center)
# correct area element size
for idx, d in enumerate(ufo.Data) :
# get normalized connection vector from phase center towards the element data point
p = np.asarray([d.v.x, d.v.y, d.v.z]) - center
R = np.sqrt(np.dot(p, p))
p = p/R
# get normalized normal vector of the element data point
n = np.asarray([d.n.x, d.n.y, d.n.z])
n = n/np.sqrt(np.dot(n,n))
# dot product gives cosine of the angle between direct connection and normal vector. Area element must have the size corresponding to a orthogonal incident angle
ufo.Data[idx].df *= np.absolute(np.dot(n, p))/ R / R * newR * newR
ufo.Data[idx].n.x = -p[0]
ufo.Data[idx].n.y = -p[1]
ufo.Data[idx].n.z = -p[2]
newPos = newR * p
ufo.Data[idx].v.x = newPos[0]
ufo.Data[idx].v.y = newPos[1]
ufo.Data[idx].v.z = newPos[2]
#####################################
# generalized functions to create a standard cassegrain or gregorian telescope with a parabolic on axis main dish
[docs]
def getAngularSizes (self, diam, distance, raster, conicData, focalPointDistance = None) :
"""
Iteratively calculate angular size of a conical section from its linear diameter. This includes the curvature of the surface.
Parameters:
diam (float): diameter [mm] of the object (resulting in an angular size in theta direction)
distance (float): distance [mm] of the origin to the conical surface vertex
raster (float): linear raster size [mm], gets converted to an angulat step-size
conicData (list) : conicData of the conical surface as calculated by UFO.UFOGaussCalc.calcConicData
focalPointDistance (float): distance of the conical surface reflection point to the focal point
Returns:
tuple (tuple): sizeDeg (float, degree), radius (float, mm, from origin to surface edge), rasterDeg (float, angular raster (float, degree), converged (boolean, flag)
"""
if focalPointDistance == None :
focalPointDistance = distance
converged = True
sizeDeg = np.arctan2(diam/2., distance) # calculate starting value for the angular size
radius = abs(distance) # set a starting radius
v, n, df, th = UFODataClass.UFOData._conicS ( sizeDeg, 0., conicData ) # calculate conical surface in cartesian coordinated from radius and angle
count = 0 # set conter to stop iterations if not converting
while abs( abs( v.x ) - diam/2. ) > 0.1 : # loop iterations until diameter is reached (degree of fredom is the radius)
count += 1 # increase iteration counter
radius = radius * 2 * abs( v.x ) / diam # adapt radius using the quotient of given and actual diameter
sizeDeg = np.arctan2(diam/2., radius) # calculate new angular size
v, n, df, th = UFODataClass.UFOData._conicS ( sizeDeg, 0., conicData ) # re-calculate conical data for the point
if count > 500 : # check non converging state
converged = False
break
sign = 1.
if conicData[2] > 1. : # check for hyperbola for the overall readius
sign = -1.
sizeDeg *= 180 / np.pi # transfer radians to degree
# radius is fictive radius for conic section, the radius for a sphere acting as spill-over region must be calculated from edge itself.
vz, nz, df, th = UFODataClass.UFOData._conicS ( 0, 0., conicData )
radius = np.sqrt( v.x**2. + (sign*np.abs(focalPointDistance) + vz.z - v.z)**2. )
rasterDeg = np.arctan2( raster, np.abs(distance) ) * 180. / np.pi # calculate angular raster size
return sizeDeg, radius, rasterDeg, converged
[docs]
def makeSecondary ( self, wavelength, SecR2, SecF, CalcSpill = True) :
"""
Create the secondary mirror out of the data set in the class.
Parameters:
wavelength (float): wavelength in mm for which the rastering of the surface is calculated. This is the lower end of the usable frequency range for simulations
SecR2 (float): distance to phase center [mm], or distance to nominal geometric focus point (this is equivalent for geometric optics)
SecF (float): focal length of the secondary mirror [mm]
CalcSpill (boolean): indicating if a spill over region is calculated and inserted into the simulation setup.
"""
# calculate rastering of the secondary mirror
SecRaster = wavelength * self.Raster / _TELESCOPE_RASTER / 2. # granularity of the surface = half wavelength
if SecRaster > self.secDiam / 300 : SecRaster = self.secDiam / 300
# Data calculated from above just to make use of the "makeConicSection" method of the UFO calss -- no direct access to the hyperbola data ...
SecWL = 5 # wavelength at which the data for the hyperbola are calculated
SecD2 = 3000. # arbitary choosen value
SecW2 = np.sqrt(np.sqrt(SecWL**2 * (SecR2*SecD2-SecD2**2)/(np.pi**2)))
# calculate missing gaussian optics data from existing
if ( SecF < 0 ) :
w1, w2, d1, d2, f = complete(SecWL, w1 = None, w2 = SecW2, d1 = None, d2 = SecD2, f = SecF)
else :
w1, w2, d1, d2, f = complete(SecWL, w2 = None, w1 = SecW2, d2 = None, d1 = SecD2, f = SecF)
#w1, w2, d1, d2, f = complete(SecWL, w1 = None, w2 = SecW2, d1 = None, d2 = SecD2, f = SecF)
# calculate conical data of the hyperbola
conicData = calcConicData ( w1, w2, d1, d2, SecWL, 0.)
# get angular size iteratively
if ( SecF < 0 ) :
fPlaneDist = SecR2
else :
fPlaneDist = SecR2 # -conicData[-1]-conicData[-2]
SecSizeDeg, radius, SecRasterDeg, converged = self.getAngularSizes(self.secDiam, SecR2, SecRaster, conicData, fPlaneDist )
# set conical section data for the sourounding sphere as spill-over area
conicSphere = (
radius, # a
radius, # b
0., # e
0., # c
radius, # p
0., # center x
0, # center z
180., # tangential-angle
False, # turn flag
radius/2., # focal-length
0., # reflection angle
radius, # r1
radius # r2
)
# defining a density function to reduce number of data points in the spill-over region
def density (x) :
if x < SecSizeDeg : return 1
if x < 2 * SecSizeDeg : return 4
return 16
# and make the conic section including spill-over area
secondary = UFODataClass.UFOData( log = self.log )
secondary._makeConic ( [SecSizeDeg, 85], SecRasterDeg, density, [conicData, conicSphere], Name = "Secondary-mirror", random = False, phiStepMult = 1.7)
# set up spill-over area
secSpill = deepcopy( secondary ) # data set of the spill-over region
DataSky = [] # Data points of the telescope active area
DataSpill = [] # Data points of the spill-over region at zenth direction (Termination at Tsky)
# select between spill-over and telescope area
zHeight = 0. # Height of the subreflector
for idx in range( 0, secondary.points ) : # loop all data points
d = secondary.Data[idx]
rx2 = d.v.x**2 # and calculate radial distance to mirror center
ry2 = d.v.y**2
r2 = rx2 + ry2
# check diameter and blockage
if ( r2 < self.secDiam**2. / 4. ) : # is data point on signal path?
DataSky.append( d ) # yes, add it to mirror data
dist = d.v.z - secondary.Header.v.z # and get maximum mirror height
if np.absolute(dist) > np.absolute(zHeight) :
zHeight = dist
else : # no, add data point to the spillover area
DataSpill.append( d )
# calculate best illumination from given edge-taper
z = SecR2 - zHeight
wz = np.sqrt( -10 * 2. * self.secDiam**2 / 4. / self.EdgeTaper / np.log( 10 ) )
# and set opening angle, and optical functions for waist size and focus
self.OpeningAngle = np.arctan2( self.secDiam / 2., z ) * 180. / np.pi
self.waist = lambda wl : np.sqrt( -np.sqrt( wz**4. / 4. - z**2. * wl**2. / ( np.pi**2. ) ) + wz**2. / 2. )
self.focus = lambda wl : self.FocalLength / 2. + np.sqrt( self.FocalLength**2. / 4. - np.pi**2 * self.waist( wl )**4. / ( wl**2. ) ) - self.FocalLength
# set secondary data
secondary.points = len( DataSky )
secondary.Data = ( UFODataClass.UFODataPoint * secondary.points )( *DataSky )
secondary.move( 0., 0., - np.abs( fPlaneDist ) )
secondary.Name = "Subreflector"
# add secondary to mirror-list
self.mirrorList.append( secondary )
self.mirrorName.append( secondary.Name )
if CalcSpill: # use spill-over area?
# set data of the secondary spill-over region
secSpill.points = len( DataSpill )
secSpill.Data = ( UFODataClass.UFODataPoint * secSpill.points )( *DataSpill )
if ( SecF < 0 ) :
secSpill.rotX( 180. )
else :
secSpill.move( 0., 0., - np.abs( fPlaneDist ) )
secSpill.Name = "Subreflector spill-over area"
self.spillList.append( secSpill )
self.spillStep.append( -1 ) # calculate spill-over from detector
self.spillTemp.append( self.Tsky ) # spill over terminated at Sky temperature
self.spillNames.append('sky')
# save telescope mirrors and spill-over regions
if self.usedFolder != "" :
secondary.writeFile (self.usedFolder + "secondary.tclass")
if CalcSpill :
secSpill.writeFile (self.usedFolder + "secondarySpillSky.tclass")
[docs]
def makePrimary ( self, wavelength, PrimF, primZpos, primRot = True, CalcSpill = True, blockage = -1., legWidth = -1., legRadius = -1., legDistance = -1., spillPrim = True , rotSupport = 0.) :
"""
Create the primary mirror out of the data set in the class.
Parameters:
wavelength (float): wavelength in mm for which the rastering of the surface is calculated. This is the lower end of the usable frequency range for simulations
PrimF (float): focal length of the secondary mirror [mm]
primZpos (float): z-position of the primary vertex [mm]
primRot (boolean): False = facing upwards, no rotation, True = facing downwards means rotating sky
CalcSpill (boolean): indicating if a spill over region is calculated and inserted into the simulation setup.
blockage (float): radius of central blockage [mm]; if negative, default value is calculated (size of the sub-reflector)
legWidth (float): leg-half-width [mm]; if negative, default value of 0. mm is taken
legRadius (float): radius at which the support legs hit the subreflector [mm]; if negative, primary radius is taken
legDistance (float): radial distance at which the legs pass the primary focal point [mm], if negative, subreflector radius is taken.
rotSupport (float): Angle of the support of the secondary (0 : horizontal and vertical)
"""
# select if spill-over is calculated from detector or from secondary mirror. If primary is rotated it can't be the detector
spillStep = -1
if primRot :
spillStep = 0
if type(self.secDiam) == type(None) :
secDiam = 0.
else:
secDiam = self.secDiam
# structure data:
if blockage < 0. :
blockage = secDiam / 2. # radius of the central blockage
if legWidth < 0. :
legWidth = 0. # leg-width (half-width)
if legRadius < 0. :
legRadius = self.Diameter/ 2. # radius at which legs hit primary
if legDistance < 0. :
legDistance = secDiam / 2. # distance leg -> focus center
# helper values for the support structure:
legHeight = legRadius**2. / (4. * PrimF )
mb = -( PrimF - legHeight ) / ( legRadius - legDistance ) # gradient of the support leg
bb = PrimF - mb * legDistance
# calculate rastering of the primary mirror
PrimRaster = wavelength * self.Raster
if PrimRaster > self.Diameter / 300 :
PrimRaster = self.Diameter / 300
if PrimRaster < self.Diameter / 2000 : # limit telescope memory consumption and calculation time
PrimRaster = self.Diameter / 2000
# set parabola data as conical data
conicPrim = (
PrimF, # a
0, # b
1., # e
0, # c
2 * PrimF, # p
0., # center x
PrimF, # center z
180., # tangential-angle
False, # turn flag
PrimF, # focal-length
0., # reflection angle
PrimF, # r1
1e99 # r2
)
# get angular size iteratively
PrimSizeDeg, radiusS, PrimRasterDeg, converged = self.getAngularSizes( self.Diameter, -PrimF, PrimRaster, conicPrim )
# set conical section data for the sourounding sphere
conicSphere = (
1.*radiusS, # a
1.*radiusS, # b
0., # e
0., # c
1.*radiusS, # p
0., # center x
PrimF, # center z
180., # tangential-angle
False, # turn flag
PrimF, # focal-length
0., # reflection angle
radiusS, # r1
radiusS # r2
)
# set density function
def density (x) :
if x < PrimSizeDeg : return 1
elif x < 75 : return 10
else : return 15
# make conic section together with spill-over sphere
primary = UFODataClass.UFOData( log = self.log )
primary._makeConic ( [ PrimSizeDeg, 90 ], PrimRasterDeg, density, [ conicPrim, conicSphere ], Name = "Primary-mirror", phiStepMult = 1.7)
DataSky = [] # Data points of the telescope active area
DataSpill = [] # Data points of the ground spill-over region at zenth direction (Termination at Tgrd)
DataSpillInner = [] # Data points of the inner blockage (Termination at Tskatter)
DataSpillSky = [] # Data points of the spill-over region dumped on sky (elevation dependent
zMin = PrimF * 2.
zMax = 0.
for idx in range(0, primary.points) : # loop all points of the parabola and check in which area they are (spill or telescope)
d = primary.Data[idx]
px = np.cos( rotSupport / 180. *np.pi ) * d.v.x - np.sin( rotSupport / 180. *np.pi ) * d.v.y
py = np.sin( rotSupport / 180. *np.pi ) * d.v.x + np.cos( rotSupport / 180. *np.pi ) * d.v.y
rx2 = px**2 # calculate radialdistance to bore-sight
ry2 = py**2
r2 = rx2 + ry2
out = False # output flag: True -- part of the active surface; False -- part of a spill-overarea
# check diameter and blockage
if ( r2 < self.Diameter**2. / 4. ) and ( r2 > blockage**2. ) :
# check inner leg part
if ( ( r2 > legRadius**2 ) or ( ( rx2 >= legWidth**2 ) and ( ry2 >= legWidth**2 ) ) ) :
if r2 < legRadius**2 :
out = True
else :
#check shadowing
m = ( d.v.z - PrimF ) / np.sqrt( r2 )
xs = ( bb - PrimF ) / ( m - mb )
ys = m * xs
w = np.sqrt( ( r2 + ( PrimF - d.v.z )**2 ) / ( xs*xs + ys*ys ) ) * legWidth
if ( rx2 > w * w ) and ( ry2 > w * w ) :
out = True
if out : # data point on signal path?
DataSky.append( d ) # add it to the appropriate data set
if d.v.z < zMin : zMin = d.v.z # and getmaximum and minimum height
if d.v.z > zMax : zMax = d.v.z
else : # data point is in the spill-area:
if ( r2 < self.Diameter**2. / 4. ) : # data point over the dish ==> scattered to sky ==> use Tscatter
DataSpillInner.append( d )
else : # data point outside the effective dish
# calculate angle towards horizon (elevation dependent)
angle = np.arctan2( d.v.z - PrimF, d.v.x ) * 180. / np.pi - self.Elevation + 90.
if angle > 180. :
angle -= 360.
# check if terminated on sky
if ( angle <= 0. ) :
DataSpill.append( d )
else:
DataSpillSky.append( d )
primary.Data = []
spill = deepcopy( primary )
spillInner = deepcopy( primary )
spillSky = deepcopy( primary )
# set telescope mirror
primary.points = len( DataSky )
primary.Data = ( UFODataClass.UFODataPoint * primary.points )( *DataSky )
# and add RMS-uncertainty
primary.surfaceError( self.RMS )
if primRot :
primary.rotX ( 180 )
primary.move ( 0., 0., primZpos )
#add telescope to simulation lists
self.mirrorList.append( primary )
self.mirrorName.append( primary.Name )
if CalcSpill :
if spillPrim :
# set ground spill region
spill.points = len( DataSpill )
spill.Data = ( UFODataClass.UFODataPoint * spill.points )( *DataSpill )
if primRot :
spill.rotX ( 180 )
spill.move ( 0., 0., primZpos )
if spill.points > 0 : #add spill-region to spill-over list
spill.Name = "Primary ground terminated spill-over area"
self.spillList.append( spill )
self.spillStep.append( spillStep )
self.spillTemp.append( _Tgrd )
self.spillNames.append( 'ground' )
# set sky spill region
spillSky.points = len( DataSpillSky )
spillSky.Data = ( UFODataClass.UFODataPoint * spillSky.points )( *DataSpillSky )
if primRot :
spillSky.rotX ( 180 )
spillSky.move ( 0., 0., primZpos )
if spillSky.points> 0 : #add spill-region to spill-over list
spillSky.Name = "Primary sky terminated spill-over area"
self.spillList.append( spillSky )
self.spillStep.append( spillStep )
self.spillTemp.append( self.Tsky )
self.spillNames.append( 'sky' )
# set scatter region
spillInner.points = len( DataSpillInner )
spillInner.Data = ( UFODataClass.UFODataPoint * spillInner.points )( *DataSpillInner )
if primRot :
spillInner.rotX ( 180 )
spillInner.move ( 0., 0., primZpos )
self.makeSpill ( spillInner, [ 0., 0., primZpos - PrimF ] )
if spillInner.points> 0 : #add spill-region to spill-over list
spillInner.Name = "Primary skatter area"
self.spillList.append( spillInner )
self.spillStep.append( spillStep )
self.spillTemp.append( _Tskatter )
self.spillNames.append( 'skatter' )
# write data files if requested
if self.usedFolder != "" :
primary.writeFile (Path(self.usedFolder, "primary.tclass"))
if CalcSpill :
spillInner.writeFile (Path(self.usedFolder,"primarySpillSkatter.tclass"))
if spillPrim :
spill.writeFile (Path(self.usedFolder, "primarySpillGnd.tclass"))
spillSky.writeFile (Path(self.usedFolder, "primarySpillSky.tclass"))
[docs]
def makeSky ( self, primZpos = 0., primRot = False, wavelength = _DEFAULT_WL, skySize = _SkySize ) :
"""
Create the sky representation for the simulation.
Parameters:
primZpos (float): z-position of the primary vertex [mm]
primRot (boolean): False = facing upwards, no rotation, True = facing downwards means rotating sky
"""
#sky = UFODataClass.UFOData( log = self.log )
#size = skySize * wavelength / self.Diameter * _SkyDistance # size of the 'sky'
#step = size / _SkyStep / skySize
#sky.makeFlat( ( size, size ), ( step, step ) )
conicSphere = (
_SkyDistance, # a
_SkyDistance, # b
0., # e
0., # c
_SkyDistance, # p
0., # center x
0, # center z
180., # tangential-angle
False, # turn flag
_SkyDistance/2., # focal-length
0., # reflection angle
_SkyDistance, # r1
_SkyDistance # r2
)
SkySizeDens = _SkySize * wavelength / self.Diameter *180. / np.pi / 2.
SkySizeDeg = skySize * wavelength / self.Diameter *180. / np.pi / 2.
SkyRasterDeg = wavelength / self.Diameter / _SkyStep * 180. / np.pi
# defining a density function to reduce number of data points in the spill-over region
def density (x) :
if x < SkySizeDens : return 1
if x < 2 * SkySizeDens : return 4
if x < 4 * SkySizeDens : return 8
if x < 10 * SkySizeDens : return 16
return 100
# and make the conic section including spill-over area
sky = UFODataClass.UFOData( log = self.log )
sky._makeConic ( [SkySizeDeg], SkyRasterDeg, density, [conicSphere], Name = "Sky", random = False, phiStepMult = 1.)
if primRot :
sky.rotX ( 180. )
sky.move( 0, 0, primZpos -_SkyDistance )
else :
sky.move( 0, 0, primZpos +_SkyDistance )
sky.Name = "Sky"
# and add it to mirrorList
self.mirrorList.append( sky )
self.mirrorName.append( sky.Name )
# save telescope mirrors and spill-over regions
if self.usedFolder != "" :
sky.writeFile ( Path(self.usedFolder, "sky.tclass") )
[docs]
def GregorianOrCass ( self,
wavelength,
CalcSpill = True,
RMS = 0.,
spillPrim = True,
TelRMS = 0.,
TelDiam = 12000.,
TelTaper = -12.,
TelPrimF = 4800.,
TelSecDiam = 750.,
TelSecR1 = -294.,
TelSecR2 = 5880.,
TelBlockage = -1,
TelLegWidth = 50.,
TelLegRadius = 3500.,
TelLegDistance = -1,
TelRotSupport = 0.
) :
"""
Creates a generalized Gregorian or Cassegrain telescope.
Parameters:
wavelength (float): wavelength to adjust gridding
CalcSpill (boolean): Shall the spill-over area be calculated and set for simulation
RMS (float): RMS of the surface (< 0 : uses telescope RMS as measured; >= 0: uses given number for the surface RMS)
spillPrim (boolean): flag if primary mirror outer spill-over is used or not (mainly because of calculation time for higher frequencies)
**kargs (???): as defined in constants for the telescope dictionaries
"""
self.printOut("Creating a Gregorian or Cassegrain telescope", self.log.FLOW_TYPE)
if RMS < 0 :
self.RMS = TelRMS # actual surface RMS of the primary
else :
self.RMS = RMS
# Mirror details
self.Diameter = TelDiam # diameter of the pimary mirror
self.EdgeTaper = TelTaper # Edge-Taper for optimum illumination (here as of the spec sheet)
# Secondary mirror
self.secDiam = TelSecDiam # diameter of the secondary mirror
TelSecMag = np.absolute( TelSecR2 / TelSecR1 ) # Magnification of the secondary
TelSecF = TelSecR1 * TelSecR2 / ( TelSecR2 + TelSecR1 ) # focal length of the secondary mirror
self.FocalLength = TelSecMag * TelPrimF # focal length of the cassegrain focus
# make secondary mirror
self.makeSecondary ( wavelength, TelSecR2, TelSecF, CalcSpill = CalcSpill )
# make primary mirror
primPos = self.mirrorList[-1].Header.v.z + TelSecR1 + TelPrimF
self.makePrimary ( wavelength, TelPrimF, primPos, primRot = True, CalcSpill = CalcSpill,
blockage = TelBlockage, legWidth = TelLegWidth, legRadius = TelLegRadius , legDistance = TelLegDistance, spillPrim = spillPrim, rotSupport = TelRotSupport)
# make sky
self.makeSky ( -TelSecR2 - TelSecR1 + TelPrimF, primRot = True )
self.printOut ("done", self.log.DEBUG_TYPE, False)
[docs]
def Parabolic ( self,
wavelength,
CalcSpill = True,
RMS = 0.,
spillPrim = True,
TelRMS = 0.,
TelDiam = 100000.,
TelTaper = -12.,
TelPrimF = 30000.,
TelBlockage = -1,
TelLegWidth = -1,
TelLegRadius = -1,
TelLegDistance = -1,
TelRotSupport = 0.
) :
"""
Creates a generalized Parabolic main dish (primary focus telecope).
Parameters:
wavelength (float): wavelength to adjust gridding
CalcSpill(boolean): Shall the spill-over area be calculated and set for simulation
RMS (float): RMS of the surface (< 0 : uses telescope RMS as measured; >= 0: uses given number for the surface RMS)
spillPrim (boolean): flag if primary mirror outer spill-over is used or not (mainly because of calculation time for higher frequencies)
**kargs (???): as defined in constants for the telescope dictionaries
"""
self.printOut("Creating a Parabolic primary focus telescope", self.log.FLOW_TYPE)
if RMS < 0 :
self.RMS = TelRMS # actual surface RMS of the primary
else :
self.RMS = RMS
# Mirror details
self.Diameter = TelDiam # diameter of the pimary mirror
self.EdgeTaper = TelTaper # Edge-Taper for optimum illumination (here as of the spec sheet)
self.secDiam = None # diameter of the secondary mirror (not existant)
self.FocalLength = TelPrimF # focal length of the telescope
# make primary mirror
self.makePrimary ( wavelength, TelPrimF, -TelPrimF, primRot = False, CalcSpill = CalcSpill,
blockage = TelBlockage, legWidth = TelLegWidth, legRadius = TelLegRadius , legDistance = TelLegDistance, spillPrim = spillPrim, rotSupport = TelRotSupport)
# calculate best illumination from given edge-taper
z = self.FocalLength - self.Diameter**2. / 16. / self.FocalLength # focal-length - parabola height at the edge
wz = np.sqrt( -10 * 2. * self.Diameter**2 / 4. / self.EdgeTaper / np.log( 10 ) )
# and set opening angle, and optical functions for waist size and focus
self.OpeningAngle = np.arctan2( z, self.Diameter / 2. )
self.OpeningAngle = 90. - 180. * self.OpeningAngle / np.pi
self.waist = lambda wl : np.sqrt( -np.sqrt( wz**4. / 4. - z**2. * wl**2. / ( np.pi**2. ) ) + wz**2. / 2. )
self.focus = lambda wl : self.FocalLength / 2. + np.sqrt( self.FocalLength**2. / 4. - np.pi**2 * self.waist( wl )**4. / ( wl**2. ) ) - self.FocalLength
# make sky
self.makeSky ( -TelPrimF, primRot = False )
self.printOut("done", self.log.DEBUG_TYPE, False)