#!/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
#############################################################################
# - clean-up code (frequently and regularly, please)
# - improve readability and add more comments (frequently and regularly)
#
# single task-changes / add-ons
#
# - write conic section in python or include it with c-type part
# - mostly done, needs more testing
#
# - add phase-corrected mirror and more functionality to shape the surfaces
# - add reflection losses for lens simulations (partly done, needs testing)
# - lenses must be newly implemented (to be done!!)
# - permittivity must be implemented correctly to handle reflection and transmitted fields
#
# - improve overall structure (still depending too much on old code)
#
#############################################################################
# import libs
#############################################################################
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.mlab as ml
import subprocess
from copy import copy, deepcopy
_GridData = False
try:
from scipy.interpolate import griddata
from scipy.optimize import curve_fit
_GridData = True
except:
_GridData = False
import LogClass
from Simulation import FarFieldDataClass
from Simulation.constants import _ETA, _Tiso, _I0iso, _Impedance, _GaussBeamDiam
from Simulation import getColors
from FileIO import FileDataClass
# import UFO-libs
from Simulation.UFO.UFO_CCode import _UFOVector, _UFOMatrix, UFODataPoint, UFOHeaderPoint, _cCalc_CurrentDistFull, _cMove, _cTransform, _cCalc_setLog
from Simulation.UFO.UFOGaussCalc import calcConicData, complete
from Simulation.UFO.UFOFileIO import _UFOFile
#############################################################################
# constants
#############################################################################
_CLASSNAME = "UFOData"
# UFO interface constants
_SURFACE = 8 # UFO surface file has 8 columns
_CURRENT = 20 # UFO current file has 20 columns
_INTENSITY = 10 # UFO intensity file has 8 columns and no header-line
#
_STEPLENS = 0.1 # relative stepwidth in iterative lens calculation (actual step-width = _STEPLENSE/f)
#############################################################################
# classes
#############################################################################
[docs]
class UFOData :
"""
View the UFO documentation page for an overview of the UFO modules.
The name UFO is unimportant. It precedes the PAF Frontend simulator by many years and was initially a standalone project.
The PAF Frontend simulator uses it for its physical optics simulation features.
This class handles UFO data (current and structure).
The data represents all the surfaces within the simulated environment of the PAF and stores their optical and electrical properties
The environment includes the reflective dishes of the telescope, the sky, spill-over regions and more.
"""
#####################################
# internal data and init method
# init method
def __init__ ( self, log = None ):
"""
Initialize class.
Parameters:
log (LogClass): Logging Class used by the calling application
"""
self.setLog( log )
# Coordinate system definitions (used to rotate and shift the data in 3D space
self.Header = UFOHeaderPoint()
self.DataID = "" # file name or other ID of the data set
self.Data = [] # data array
self.points = 0 # size of the data array
self.Type = _SURFACE; # file type (implies number of columns)
self.Wavelength = None # if simulation results are present, the wavelength is stored here (is not saved)
self.I0 = 0. # if set, it can be used to determine power of the final beam
self.er = 1. # if set, this permittivity is used for simulations from this object to the next
self.tanD = 0. # if set, this loss-factor is used for simulations towards this object
self.SurfaceType = "Unknown" # if set, it tells the surface type represented by the point-cloud (e.g. "Aperture", indicting result must be inverted)
self.FocalLength = None # if set, this represents the focal length of the object
self.RefAngle = 0. # if set, this represents the reflection angle of the object
self.Name = "--" # if set, this represents the name of the object
[docs]
def setLog ( self, newlog = None ) :
"""
Set new log class
"""
if newlog == None:
self.log = LogClass.LogClass( _CLASSNAME )
else:
self.log = deepcopy( newlog ) # logging class used for all in and output messages
self.log.addLabel( self )
[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)
#####################################
# File operations
[docs]
def readFile ( self, filename ):
"""
Read UFO file.
Parameters:
filename (string): filename of the UFO file (structure or current)
"""
fileNameData = FileDataClass.FileNameData(filename) # parse filename
self.DataID = fileNameData.baseName # and extract base-name as ID
FD = FileDataClass.FileData( headerProfile = _UFOFile, log = self.log ) # initialize file data with UFO profile
Header, self.Data = FD.readData( filename ) # read file header and data
self.points = len(self.Data) # get number of data-points
self.Data = (UFODataPoint*self.points)(*self.Data) # and translate all data points to memory array
# extract header point
self.Header = UFOHeaderPoint() # extract center location and orientation
self.Header.v = _UFOVector( Header['CENTER:'][0], Header['CENTER:'][1], Header['CENTER:'][2] )
self.Header.o = _UFOVector( Header['ORIENTATION_X:'][0], Header['ORIENTATION_X:'][1], Header['ORIENTATION_X:'][2] )
self.Header.n = _UFOVector( Header['ORIENTATION_Z:'][0], Header['ORIENTATION_Z:'][1], Header['ORIENTATION_Z:'][2] )
# set header data # set header data from header dictionary
self.Type = Header['FILESTATE:'][0]
self.Wavelength = Header['WAVELENGTH:'][0]
if self.Wavelength == 'None' :
self.Wavelength = None
else :
self.Wavelength = float( self.Wavelength )
self.I0 = Header['I0:'][0]
self.er = Header['ER:'][0]
self.tanD = Header['TAND:'][0]
self.SurfaceType = Header['TYPE:'][0]
self.FocalLength = Header['F:'][0]
if self.FocalLength == 'None' :
self.FocalLength = None
else :
self.FocalLength = float( self.FocalLength )
self.RefAngle = Header['REFANGLE:'][0]
self.Name = Header['NAME:'][0]
[docs]
def writeFile ( self, filename ) :
"""
Write UFO file.
Parameters:
filename (string): filename of the UFO file (structure or current, depending on the internal Type)
"""
# set header dictionary
Header = {
'WAVELENGTH:' : [ str(self.Wavelength) ],
'I0:' : [ self.I0 ],
'ER:' : [ self.er ],
'TAND:' : [ self.tanD ],
'TYPE:' : [ self.SurfaceType ],
'F:' : [ str(self.FocalLength) ],
'REFANGLE:' : [ self.RefAngle ],
'NAME:' : [ self.Name ],
'FILESTATE:' : [ self.Type ],
'CENTER:' : [ self.Header.v.x, self.Header.v.y, self.Header.v.z ],
'ORIENTATION_X:' : [ self.Header.o.x, self.Header.o.y, self.Header.o.z ],
'ORIENTATION_Z:' : [ self.Header.n.x, self.Header.n.y, self.Header.n.z ],
}
FD = FileDataClass.FileData( headerProfile = _UFOFile, log = self.log ) # initialize file data with UFO profile
FD.writeData( filename, data = self.Data, header = Header )
[docs]
def getData ( self, onlyPos = False ) :
"""
Returns the internal data as numpy arrays.
Parameters:
onlyPos (boolean): if True, only the position array is returned
Returns:
tuple (tuple): of Numpy array(s)
Numpy array(s):
1) position data of all mirror points
2) area element sizes (only if onlyPos == False)
3) complex E-field (only if onlyPos == False)
4) complex current distribution (only if onlyPos == False)
"""
self.printOut("Returning data as numpy arrays :", self.log.FLOW_TYPE)
E_ = np.zeros( (self.points, 3 ), dtype = np.complex128)
I_ = np.zeros( self.points, dtype = np.float64 )
v_ = np.zeros( (self.points, 3 ), dtype = np.float64 )
df_ = np.zeros( self.points, dtype = np.float64 )
for idx in range(0, self.points) :
v_ [idx] = np.array([self.Data[idx].v.x, self.Data[idx].v.y, self.Data[idx].v.z])
df_[idx] = self.Data[idx].df
if self.Type == _SURFACE :
E_ = np.array( [0. + 0 * 1j, 0. + 0 * 1j, 0. + 0 * 1j] )
I_ = 0.
else:
E_ [idx] = np.array([self.Data[idx].Er.x+1j*self.Data[idx].Ei.x, self.Data[idx].Er.y+1j*self.Data[idx].Ei.y, self.Data[idx].Er.z+1j*self.Data[idx].Ei.z])
I_ [idx] = np.absolute(self.Data[idx].I)
if onlyPos :
return v_
return v_, E_, I_, df_
#####################################
# Field import/export or calculation
[docs]
def makeFarField ( self, wavelength = 100., cx = 0., cy = 0., cz = 0. ) :
"""
Does not work properly!!!
Tries to create a far-field data class file out of the UFO data present.
(cx, cy, cz: best center of a sphere, but all other points work as well. Wrong data can be corrected with FarFieldDataClass in post processing.)
Parameters:
wavelength (float): required if phase should be retrieved in an absolute manner
cx (float): Coordinate of phase center of the farfield data.
cy (float): Coordinate of phase center of the farfield data.
cz (float): Coordinate of phase center of the farfield data.
"""
self.printOut("Creating Far-Field data from UFO current file", self.log.FLOW_TYPE)
if self.Type != _CURRENT :
self.printOut("no field data available ( instance is not of type 'CURRENT'", self.log.ERROR_TYPE, False)
self.printOut("using :", self.log.DATA_TYPE, False)
self.printOut("wavelength : %.2fmm" % wavelength, self.log.DATA_TYPE, False)
self.printOut("center : (%.2fmm, %.2fmm, %.2fmm)" % (cx, cy, cz), self.log.DATA_TYPE, False)
self.printOut("extract field data", self.log.USER_TYPE, False)
ff = FarFieldDataClass.FarFieldData(wavelength, log = self.log ) # establish far-field data
center = np.asarray([cx, cy, cz]) # get center-vector
zeroPos = np.array([self.Header.v.x, self.Header.v.y, self.Header.v.z])
refDist = np.sqrt(np.dot(zeroPos-center, zeroPos-center)) # reference distance is distance to UFO center
# creating interpolation-function
x = np.zeros( self.points )
y = np.zeros( self.points )
E = np.zeros( ( self.points, 3 ), dtype = complex )
maxR = 0.
for idx in range(0, self.points) :
# caclulate pointing-vector
r = np.array( [ self.Data[idx].v.x, self.Data[idx].v.y, self.Data[idx].v.z ] ) - center
x[idx] = r[0]
y[idx] = r[1]
# calculate phi and theta
Phi = np.arctan2( r[1], r[0] )
if Phi < 0. : Phi += 2.*np.pi
Theta = np.arctan2( np.sqrt(r[0]**2.+r[1]**2.), r[2] )
# scale intensity of the E-field and match phase to far-field reference distance
E_ = np.array( [ self.Data[idx].Er.x+1j*self.Data[idx].Ei.x, self.Data[idx].Er.y+1j*self.Data[idx].Ei.y, self.Data[idx].Er.z+1j*self.Data[idx].Ei.z ] )
E[idx] = E_ / refDist * ff.ReferenceDist * np.exp( -1j*np.pi*2.*(refDist - ff.ReferenceDist)/ff.Wavelength )
# get maximum distance to center to determine theta-range
R = r[0]**2+r[1]**2
if R > maxR: maxR = R
maxR = np.sqrt(maxR)
# create far-field data array
self.printOut("create far-field data set", self.log.DEBUG_TYPE, False)
ff.FarField = []
# get an approximation of the median point-difference
diff = x**2. + y**2
idxList = np.argsort(diff)
unit = 2.*diff[idxList[-1]]/np.sqrt(len(x))
phiRange = np.linspace(0, 360, 361)
thetaRange = np.linspace(0, np.arctan2(maxR, refDist)*180./np.pi, 180)
for phi in phiRange :
line = []
for theta in thetaRange :
r = np.tan(theta/180*np.pi)*refDist
xpos = r*np.cos(phi/180.*np.pi)
ypos = r*np.sin(phi/180.*np.pi)
if (xpos >= np.min(x)) and (xpos <= np.max(x)) and (ypos >= np.min(y)) and (ypos <= np.max(y)) :
# calculate differences to available points
diff = (x-xpos)**2. + (y-ypos)**2
idxList = np.argsort(diff)
# check if it is close to one point
if diff[idxList[0]] < unit/100. :
EField = E[idxList[0]] # then use this point
else :
EField = np.zeros(3, dtype = complex) # otherwise interpolate between next 4 neighbors
norm = 0.
for i in range(0, 4):
EField += E[idxList[i]]/diff[idxList[i]]
norm += 1./diff[idxList[i]]
EField = EField/norm
ST = np.asarray([
[np.sin(theta)*np.cos(phi), np.sin(theta)*np.sin(phi), np.cos(theta)],
[np.cos(theta)*np.cos(phi), np.cos(theta)*np.sin(phi), -np.sin(theta)],
[-np.sin(phi), np.cos(phi), 0. ]
]).T
Ep = np.dot(ST, EField.conj().T)
line.append([phi, theta, Ep[2], Ep[1]])
ff.FarField.append(line)
return ff
[docs]
def makeWaist ( self, waist, polAngle = 0., stepPerWaist = 25, size = [], wavelength = None, new = False) :
"""
"""
if (len(self.Data) == 0 ) or ( new ) :
if isinstance(size, float) :
size = [ size, size ]
if len(size) < 1 :
size = [waist*5, waist*5.]
self.printOut("Creating a waist of %.2fmm on a flat mirror" % waist, self.log.FLOW_TYPE)
self.makeFlat( dim = size, step = [waist/stepPerWaist, waist/stepPerWaist], Name = "Waist (%.2f mm)" % waist )
else :
self.printOut("Creating a waist of %.2fmm on an existing mirror" % waist, self.log.FLOW_TYPE)
self.Type = _CURRENT
center = np.array( [ self.Header.v.x, self.Header.v.y, self.Header.v.z ] )
ax, ay, az = self.getRot()
ax = -ax/180*np.pi
ay = -ay/180*np.pi
az = -az/180*np.pi
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)]
])
rz = np.asarray([
[np.cos(az), -np.sin(az), 0.],
[np.sin(az), np.cos(az), 0.],
[ 0., 0., 1.]
])
matrix = np.dot(rx, np.dot( ry, rz ) )
for idx in range( 0, self.points ) :
r = np.array( [ self.Data[idx].v.x, self.Data[idx].v.y, self.Data[idx].v.z ] ) - center
r = np.sqrt( np.dot( r, r ) )
E = 1000. * np.sqrt( np.exp( -2. * r**2. / ( waist**2. ) ) / ( 0.5 * np.pi * waist**2. ) )
Er = E*_ETA/(4*np.pi) * np.dot(matrix, np.asarray( [np.cos(np.pi*polAngle/180.), np.sin(np.pi*polAngle/180.), 0.] ) )
Jr = -2 / _ETA * Er
self.Data[idx].Er = _UFOVector( Er[0], Er[1], Er[2] )
self.Data[idx].Jr = _UFOVector( Jr[0], Jr[1], Jr[2] )
self.Data[idx].Ei = _UFOVector(0., 0., 0.)
self.Data[idx].Ji = _UFOVector(0., 0., 0.)
self.Data[idx].I = (0.5/_ETA) * (self.Data[idx].Er.x**2 + self.Data[idx].Er.y**2 + self.Data[idx].Er.z**2)
self.Wavelength = wavelength
self.setI0(_I0iso, setOnly = False) # calculate excitation current
self.SurfaceType = "Horn"
self.FocalLength = None
self.RefAngle = 0.
[docs]
def makeAperture ( self, diam, polAngle = 0., stepPerDiam = 100, tiltAngle = 0., Name = "") :
"""
Creates an apperture by using a flat surface which, afer calculating surface currents, is turned by 180 degree.
Hence the apperture is a mirror which will be turned to radiate into the original light direction.
"""
self.printOut("Creating an element representing an apperture", self.log.FLOW_TYPE)
self.makeFlat( dim = [diam, diam], step = [diam/stepPerDiam, diam/stepPerDiam], Name = Name )
self.rotY( tiltAngle )
Data = []
center = np.array([self.Header.v.x, self.Header.v.y, self.Header.v.z])
for idx, d in enumerate( self.Data ) :
r = np.sqrt( d.v.x**2. + d.v.y**2. )
if r <= diam/2. :
d.df = d.df * np.cos( tiltAngle / 180. * np.pi )
d.n.x = 0.
d.n.y = 0.
d.n.z = 1.
Data.append( self.Data[idx] )
self.points = len( Data )
self.Data = ( UFODataPoint * self.points )( *Data )
self.printOut(" ending with in total %d raster points." % self.points, self.log.DEBUG_TYPE)
self.SurfaceType = "Aperture"
self.FocalLength = None
self.Header.n = _UFOVector(0., 0., -1.)
self.Header.o = _UFOVector(1., 0., 0.)
self.RefAngle = 180.
#####################################
# create mirror or lens
[docs]
def makeParabola ( self, f, dim =[10., 10.], step = [1., 1.], height = -1, Name = "" ) :
"""
Creates a parabolic mirror (roration symmetric) with opening in positive z-direction).
Parameters:
f (float): optical focal-length (f=1/p = 1/(4*a) )
dim (list): dimensions in x-y plane (list with two entries)
step (list): stepping of the mirror
height (float): maximum hight (< 0: no limit)
Name (string): identifier string of the object (optional)
"""
self.printOut("Creating a parabolic mirror in x-y plane", self.log.FLOW_TYPE)
self.printOut("Mechanical data :", self.log.DATA_TYPE, False)
self.printOut("Dimensions (x,y) : %.2fmm x %.2fmm" % (dim[0], dim[1]), self.log.DATA_TYPE, False)
self.printOut("Rastering (x,y) : (%.2fmm, %.2fmm)" % (step[0], step[1]), self.log.DATA_TYPE, False)
self.printOut("Optical data :", self.log.DATA_TYPE, False)
self.printOut("focal-length : %.2fmm" % f, self.log.DATA_TYPE, False)
# clear data and set header
self.Data = []
self.Header = UFOHeaderPoint()
self.Header.v = _UFOVector(0., 0., 0.)
self.Header.n = _UFOVector(0., 0., 1.)
# make sure (0, 0) is in the point cloud to symmetrise the problem
startPosX = int(0.5*dim[0]/step[0])*step[0]
startPosY = int(0.5*dim[0]/step[1])*step[1]
self.Type = _SURFACE
# loop all positions
for y in np.arange(-startPosY, startPosY +1, step[0]) :
for x in np.arange(-startPosX, startPosX + 1, step[1]) :
d = UFODataPoint()
d.v = _UFOVector(x, y, 0.25*(x*x+y*y)/f)
n = np.asarray([-0.5*x/f, -0.5*y/f, 1.])
nNorm = np.sqrt(np.dot(n, n))
d.df = step[0]*step[1]*nNorm
n = n/nNorm
d.n = _UFOVector(n[0], n[1], n[2])
if (height < 0) or (d.v.z < height) :
self.Data.append(d)
self.points = len(self.Data)
self.Data = (UFODataPoint*self.points)(*self.Data)
self.SurfaceType = "Mirror"
self.FocalLength = f
self.RefAngle = 0.
self.Name = Name
self.Header.o = _UFOVector(self.Header.n.z, self.Header.n.y, -self.Header.x)
[docs]
def makeConicSection ( self, d1, d2, w1, w2, fLength, wavelength, refAngle = 0, diameter = -1., step = -1., Name = "" ) :
"""
Creates a mirror from rotated conic sections (hyperbola, parabola or ellipsoidal) from parameters according Gaussian optics.
Parameters:
d1 (float): distance to Waist 1
d2 (float): distance to Waist 2
w1 (float): waist-size of waist 1
w2 (float): waist-size of waist 2
fLength (float): focal-length of the mirror
wavelength (float): wavelength being used for the calculations
refAngle (float): reflection angle of the mirror
diameter (float): diameter of the projected mirror (is automatically calculated to beam-diameter size from constants if not given)[mm]
step (float): stepping of the mirror (is automatically calculated to 0.5 lambda if not given)
Name (string): identifier string of the object (optional)
"""
# get data of the conic section from Gaussian beam parameters
conicData = calcConicData ( w1, w2, d1, d2, wavelength, refAngle)
cx = conicData[5]
cz = conicData[6]
radius = np.sqrt ( cz**2 + cx**2 )
# calculate step if missing
if step < 0 :
step = wavelength / 2.
# calculate diameter from _GaussBeamDiam
if diameter < 0 :
zr = np.pi * w1**2 / wavelength
wz = w1 * np.sqrt( 1 + ( d1 / zr )**2 )
diameter = _GaussBeamDiam * wz
thMax = np.arctan2( diameter / 2, radius ) *180. / np.pi
thStep = np.arctan2( step, radius ) *180. / np.pi
self.printOut("Creating mirror out of a rotated conical section using Gaussian beam parameters", self.log.FLOW_TYPE)
self.printOut("Gaussian optics parameters:", self.log.DATA_TYPE, False)
self.printOut("Wavelength : %.2fmm" % wavelength, self.log.DATA_TYPE, False)
self.printOut("Input : d1 = %.2fmm, w1 = %.2fmm, r1 = %.2f mm" % ( d1, w1, conicData[11] ), self.log.DATA_TYPE, False)
self.printOut("Output : d2 = %.2fmm, w2 = %.2fmm, r2 = %.2f mm" % ( d2, w2, conicData[12] ), self.log.DATA_TYPE, False)
self.printOut("focal-length : %.2f mm" % fLength, self.log.DATA_TYPE, False)
self.printOut("Reflection angle : %.2f deg" % refAngle, self.log.DATA_TYPE, False)
self.printOut("Mechanical data :", self.log.DATA_TYPE, False)
self.printOut("Diameter : %.2fmm ==> %.2f deg" % ( diameter, thMax * 2.), self.log.DATA_TYPE, False)
self.printOut("Rastering : %.2fmm ==> %.2f deg" % ( step, thStep ), self.log.DATA_TYPE, False)
self._makeConic ( [thMax], thStep, lambda x : 1., [conicData], Name = Name )
@staticmethod
def _conicS ( theta, phi, conicData ) :
"""
Calculate a point of of the canonical section including normal vector.
Parameters:
theta (float): polar angles
phi (float): polar angles
conicData (list): list with conicData
Returns:
tuple (tuple): position vector as _UFOVector, normal vector as _UFOVector, radius squared as float, th (Theta)
"""
# extract data fromconic-data set
a = conicData[0] # major half-axis
b = conicData[1] # minor half-axis
e = conicData[2] # excentricity
c = conicData[3] # conic constant
p = conicData[4] # focal parameter
cx = conicData[5] # 2d mirror center point (x-coordinate)
cz = conicData[6] # 2d mirror center point (z-coordinate)
tanD = conicData[7]/180.*np.pi #- 0.5*np.pi # angle of the mirror surface in respect to nominal conical section coordinate system
turn = conicData[8] # bool-flag, if mirror should be turned by 180deg (upside down)
f = conicData[9] # focal length
refA = conicData[10] # reflection angle
r1 = conicData[11] # distance mirror center to focal point 1
r2 = conicData[12] # distance mirror center to focal point 2
tan = -np.arctan2(cx, cz) # theta angle of the mirror center in polar coordinates of the conic section coordinate system
signRadius = 1.
if ( np.abs(r1) > np.abs(r2) ) and ( r1*r2 < 0. ) :
signRadius = -1.
tanD -= np.pi
# implement roration matrice into conic-section coordinate system
RyI = np.asarray( [
[ np.cos(tan), 0, -np.sin(tan) ],
[ 0 , 1, 0 ],
[ np.sin(tan), 0, np.cos(tan) ]
] )
# implement rotation matrice into mirror coordinate system
Ry = np.asarray( [
[ np.cos(-tanD), 0, -np.sin(-tanD) ],
[ 0 , 1, 0 ],
[ np.sin(-tanD), 0, np.cos(-tanD) ]
] )
# transform polar coordinated to cartesian coordinates
v = np.asarray( [ np.sin( theta ) * np.cos ( phi ), np.sin( theta ) * np.sin ( phi ), np.cos( theta ) ] )
# rotate to mirror system
v = np.dot( RyI, v )
# calculate local polar coordinates
th = np.arctan2( np.sqrt( v[0]**2 + v[1]**2 ), v[2] )
ph = np.arctan2( v[1], v[0] )
# get conic section radius
r = p / ( 1 + signRadius * e * np.cos( th ) )
# establish cartesian coordinates
v = np.asarray( [ r * np.sin( th ) * np.cos ( ph ) , r * np.sin( th ) * np.sin ( ph ), r * np.cos( th ) ] )
# calculate normal vector
if e == 0. : # for a sphere
R = np.sqrt( v[0]**2. + v[1]**2 + v[2]**2 )
n = np.asarray( [ -v[0] / R, -v[1] / R, -v[2] / R ] )
elif e == 1 : # for a parabola
n = np.asarray( [ -v[0] / p , -v[1] / p, -1. ] )
else : # for any other conic section
# two possiblities for the section to cut out. Depends on r1 and r2 length
# sign of the normal vector must change here
sign = 1.
if v[2] + c < 0 :
sign = -1.
# hyperbola ? If so adapt normal vector to changed
signHyper = 1.
if r1*r2 < 0 :
signHyper = -1.
sign = 1.
R = np.sqrt( v[0]**2. + v[1]**2 ) # vector length in x-y plane
if R == 0.:
n = sign * np.asarray( [ 0., 0., 1. ] )
else:
# dz = ( -R * a / ( b**2 * np.sqrt( 1 - R**2 / b**2 ) ) ).real # deviation of the 2d conic section (since b can be purely imaginary only, make result real)
# n = -1 * sign * np.asarray( [ signHyper * sign * dz * v[0] / R, signHyper * sign * dz * v[1] / R, -1.] )
if R**2 / b**2 < .5 :
dz = ( -R * a / ( b**2 * np.sqrt( 1 - R**2 / b**2 ) ) ).real # deviation of the 2d conic section (since b can be purely imaginary only, make result real)
n = -1 * sign * np.asarray( [ signHyper * sign * dz * v[0] / R, signHyper * sign * dz * v[1] / R, -1.] )
else:
dr = 1 * ( -np.abs(v[2]+c) * b / ( a**2 * np.sqrt( 1 - np.abs(v[2]+c)**2 / a**2 ) ) ).real
n = 1 * sign * np.asarray( [ signHyper * sign * v[0] / R, signHyper * sign * v[1] / R, -dr] )
# outside ellipse, r1 < 0, r2 < 0 : untested!!!!
if ( r1 < 0. ) and ( r2 < 0 ) :
n = -n
v[0] = -v[0]
n[0] = -n[0]
# normalize normal vector
n = n / np.sqrt( np.dot(n, n) )
# get normalized position vector
if np.dot( v, v ) != 0. :
vn = v / np.sqrt( np.dot( v, v ) )
else:
vn = n
# calculate area element projection including radius square ( sin( Theta), dTheta and dPhi will be multiplied in later )
df = r**2. / np.absolute( np.dot( n, vn ) )
# transform back to mirror coordinate system
# special case hyperbola: if cz is negative, use distance from other focal point to get center to zero
if cz < 0 :
cz = - signRadius * ( cz + 2. * c )
v -= np.asarray( [signRadius * cx, 0, signRadius * cz] )
v = np.dot( Ry , v )
n = np.dot( Ry , n )
if ( p > 0. ) and ( e > 0 ) and ( e != 1. ) : n *= -1.
return _UFOVector(v[0], v[1], v[2]), _UFOVector(n[0], n[1], n[2]), df, th
def _makeConic ( self, matchingPoints, stepSize, densityFunc, conicData, Name = "", random = False, phiStepMult = 1., totalAngle = 2*np.pi ) :
"""
Creates a mirror from rotated conic sections (hyperbola, parabola or ellipsoidal) from parameters according Gaussian optics.
Parameters:
matchingPoints (list): of floats, giving the angular sizes which must be met accurately. Last value gives the angular radius of the mirror [deg]
stepSize (float): nominal step-size [deg] to raster the mirror, giving the step size in theta direction. For phi-direction see also phiStepMult.
densityFunction (function): f(angle) returning float, giving a density value to modify the step size. step size used is f(theta)*stepSize
conicData (list): of lists holding the output of internal function calcConicData (must have the same length as matchingPoints
Name (string): identifier string of the object (optional)
random (boolean): randomize data points
phiStepMult (float): multiplication factor giving relative angular step-size difference between theta- and phi-steps. Phi-step = stepSize * phiStepMult
totalAngle (float): unit radiand, total phi angle used for the conic section
"""
# internal function
def dFunc ( th, densityFunc = densityFunc, random = random ) :
"""
Modify density function with random numbers.
"""
rVal = 1.
if random :
rVal = np.random.rand()*2.
return densityFunc(th)*rVal
def get_thetaRange(matchingPoints, stepSize, dFunc, random = False) :
"""
Calculate theta points in the designated range.
Parameters:
matchingPoints (list): of floats; giving the angular sizes which must be met accurately. Last value gives the angular radius of the mirror [deg]
stepSize (float): nominal step-size [deg] to raster the mirror
densityFunction (function): f(angle) returning float, giving a density value to modify the step size. step size used is f(theta)*stepSize
random (boolean): indicating if step size is randomized
"""
# main function
th = 0. # actual theta angle for appending
tRange = [] # list, (range of theta angles to be used, belonging step-size, belonging index of conicData)
s = stepSize * dFunc(th, random = False) # actual step-size
for idx, maxTh in enumerate( matchingPoints ) : # loop all marks for matching radius
while th + s/ 2. < maxTh : # check if actual theta + half step-size is still in the limit of the actual matching radius
if th == 0. : # distinguish for th = 0 because of the step-size and append theta point to range-array
tRange.append( [ th, s/2., idx ] )
else :
tRange.append( [ th, s, idx ] )
th = th + s/2. # add half step-size, coming to the edge of the actual radius
s = stepSize * dFunc(th) # get new step-size from density function
th = th + s/2. # add half of the step-size to come to next theta coordinate
# check if close to matching Edge (10% of step-size to avoide to small steps)
sCheck = (maxTh - th + s/2.) / 2. # get difference from nominal matching radius to last theta-position and divide it by 2
if sCheck > stepSize * dFunc(th, random = False) / 20. : # check resulting step-size for being larger than 10% of the actual step size at theta position
th = th - s/2. + sCheck # if so, calculate new thepa position
tRange.append( [ th, sCheck*2, idx ] ) # and append it together with the belonging step-size
th = th + sCheck # add half step-size, coming to the edge of the actual radius
s = stepSize * dFunc(th) # get new step-size from density function
th = th + s/2. # add half of the step-size to come to next theta coordinate
# transform the resulting theta-angles to radian
tRange = np.asarray(tRange)
tRange[:, 0:2] = tRange[:, 0:2] / 180. *np.pi
return tRange
# main method
self.printOut("Creating mirror composed out of %d different rotated conical sections" % len(conicData), self.log.FLOW_TYPE)
self.printOut("using matching point list [deg]: %s" %(matchingPoints), self.log.DATA_TYPE, False)
self.printOut("using step-size of : %.4f deg" %(stepSize), self.log.DATA_TYPE, False)
# chaeck data consistency
if len(matchingPoints) > len(conicData) :
self.printOut("number of conicData sets does not match number of matchingPoints, exiting", self.log.ERROR_TYPE, False)
return
tRange = get_thetaRange(matchingPoints, stepSize, dFunc, random ) # get theta values
# clear data-array
self.Data = []
dth = 0.
self.printOut("looping all theta-values (%d in total):" % ( len(tRange) ), self.log.DEBUG_TYPE, False)
for idx, thList in enumerate( tRange ) : # loop all theta values
stepPh = np.pi * dFunc( thList[0] * 180. / np.pi ) * stepSize * phiStepMult / 180. # get actual step size belonging to this theta
stepNo = int(totalAngle * np.sin(thList[0]) / stepPh + .5) # calculate number of phi-steps required
stepPhOff = ( idx % 5 - 2.5 ) * stepPh
self.printOut("theta = %.2f deg, (theta-step %d, %d phi-steps to go)" %(thList[0]*180./np.pi, idx, stepNo), self.log.FLOW_TYPE, False, False)
if stepNo == 0 : # if no phi steps are required (central point), calculate center and set header
d = UFODataPoint() # create new data point
d.v, d.n, dF, thI = self._conicS ( thList[0], 0, conicData[ int( thList[2] ) ] ) # get data point from function
d.df = np.pi * np.sin( thList[1] )**2 * dF # and adapt area size
self.Data.append(d) # append data list
self.Header = UFOHeaderPoint()
self.Header.v = deepcopy( d.v ) # and create header entry
self.Header.n = deepcopy( d.n )
self.Header.o = _UFOVector(self.Header.n.z, self.Header.n.y, -self.Header.n.x)
else : # else calculate all phi-points
stepP = 2. * np.pi / stepNo # get step size in phi direction
phi = 0. # set phi-counter to zero
for phiIdx in range(0, stepNo) : # loop all phy steps
d = UFODataPoint() # create new data point
phi = stepPhOff + stepP * phiIdx # calculate phi
if random :
dth = np.random.rand()*thList[1] / 2. - thList[1]/4
d.v, d.n, dF, thI = self._conicS ( thList[0] + dth, phi, conicData[ int( thList[2] ) ] ) # get data point from function
d.df = np.absolute( np.sin(thList[0])*dF*stepP*thList[1] ) # and adapt area size
self.Data.append(d) # append data list
self.printOut("", self.log.FLOW_TYPE, False)
self.points = len(self.Data) # set number of data points
self.Data = (UFODataPoint*self.points)(*self.Data) # transform data array to array of UFO data points for feeding the c-code
self.SurfaceType = "Mirror" # set surface type,
self.FocalLength = conicData[0][9] # focal length,
self.RefAngle = conicData[0][10] # reflection angle,
self.Name = Name # and name
self.resetPosition()
[docs]
def makeFlat ( self, dim =[10., 10.], step = [1., 1.], refAngle = 0., radius = None, Name = "" ) :
"""
Creates a flat mirror in x-y plane with normal vectors towards positive z.
Parameters:
dim (list): dimensions in x-y plane (list with two entries)
step (list): stepping of the mirror
refAngle (float): intended reflection angle at which the mirror is used (float, optional)
radius (boolean): if given, the flat mirror will have a round shape with this radius
Name (string): identifier string of the object (optional)
"""
self.printOut("Creating flat mirror in x-y plane", self.log.FLOW_TYPE)
self.printOut("Mechanical data :", self.log.DATA_TYPE, False)
self.printOut("Dimensions (x,y) : %.2fmm x %.2fmm" % (dim[0], dim[1]), self.log.DATA_TYPE, False)
self.printOut("Rastering (x,y) : (%.2fmm, %.2fmm)" % (step[0], step[1]), self.log.DATA_TYPE, False)
# clear data and set header
self.Data = []
self.Header = UFOHeaderPoint()
self.Header.v = _UFOVector(0., 0., 0.)
self.Header.n = _UFOVector(0., 0., 1.)
self.Header.o = _UFOVector(1., 0., 0.)
maxRad = 1.1*np.sqrt(dim[0]**2 + dim[1]**2)
if isinstance(radius, float) or isinstance(radius, int) :
maxRad = radius
# make sure (0, 0) is in the point cloud to symmetrise the problem
startPosX = int(0.5*dim[0]/step[0])*step[0]
startPosY = int(0.5*dim[1]/step[1])*step[1]
self.Header.columns = 0
self.Type = _SURFACE
for y in np.arange(-startPosY, startPosY +1, step[0]) :
for x in np.arange(-startPosX, startPosX + 1, step[1]) :
if x*x+y*y < maxRad**2. :
d = UFODataPoint()
d.v = _UFOVector(x, y, 0.)
d.n = _UFOVector(0., 0., 1.)
d.df = step[0]*step[1]
self.Data.append(d)
self.points = len(self.Data)
self.Data = (UFODataPoint * self.points)(*self.Data)
self.SurfaceType = "Flat"
self.FocalLength = None
self.RefAngle = refAngle
self.Name = Name
[docs]
def makeSphere ( self, radius = 1000., thetaMax = 90., stepFunction = lambda x : x + .1/60. + x/60., Name = "" ) :
"""
Creates a sphere with possible variable theta spacing.
Parameters:
radius (float): radius of the sphere
thetaMax (float): size of the sphere forward angle (0: only positive z, 90: full forward hemispher, 180: 4PI)
stepFunction (function): function to go from actual theta to next theta value default, step size increasing with theta by 1 arcmin per degree
Name (string): identifier string of the object (optional)
"""
self.printOut("Creating a sphere", self.log.FLOW_TYPE)
self.printOut("Mechanical data :", self.log.DATA_TYPE, False)
self.printOut("Radius : %.2fmm" % radius, self.log.DATA_TYPE, False)
self.printOut("maximum Theta : %.2fdeg" % thetaMax, self.log.DATA_TYPE, False)
# clear data and set header
self.Data = []
# get starting values
dtheta = stepFunction(0)
theta = dtheta
# and calculate theta = 0 point
d = UFODataPoint()
d.n = _UFOVector(0., 0., -1.)
d.v = _UFOVector(0., 0., radius)
d.df = ( radius * dtheta / 180. * np.pi )**2.
self.Data.append(d)
# create header
self.Header = UFOHeaderPoint()
self.Header.v = d.v
self.Header.n = d.n
self.Header.o = _UFOVector( 0., 0., 1. )
self.Type = _SURFACE
# loop all other points
while theta < thetaMax : # loop theta range
dphi = dtheta/np.sin(theta/180.*np.pi) # set dphi such that area element size stays nearly constant
stepsPhi = int(360./dphi)
if stepsPhi < 1 : stepsPhi = 1
dphi = 360./stepsPhi
phi = 0.
while phi < 360. : # loop phi range
d = UFODataPoint() # and calculate position, normal vector, and area size
d.n = _UFOVector(
np.cos(phi/180*np.pi) * np.sin(theta/180.*np.pi),
np.sin(phi/180*np.pi) * np.sin(theta/180.*np.pi),
-np.cos(theta/180.*np.pi)
)
d.v = _UFOVector(d.n.x*radius, d.n.y*radius, d.n.z*radius)
d.v.z *= -1.
d.df = radius*radius*dtheta*dphi*np.sin(theta/180.*np.pi)*np.pi*np.pi/180./180.
self.Data.append(d)
phi+= dphi
dthetaNew = stepFunction(theta)-theta # get next theta position step-size
theta += (dtheta+dthetaNew)/2. # next position = 0.5 x old-step size + 0.5 x new step size
# ==> data point has area belonging to actual step size
# and data point is in the center of the area
dtheta = dthetaNew
self.points = len(self.Data)
self.Data = (UFODataPoint*self.points)(*self.Data)
self.SurfaceType = "Sphere"
self.FocalLength = None
self.RefAngle = 0.
self.Name = Name
[docs]
def fromPoints ( self, fileName, inner = True, unit = 1, center = [], refAngle = 0., Name = "", Raster = 1 ):
"""
Loads a point cloud and calculates the normal vectors and area sizes.
Parameters:
fileName (string): Filename of the point cloud data file
inner (boolean): flag for the direction (inner part of the convex hull or outer)
unit (int): multiplication factor to get from file units to mm (default mm, 1000 for meter to mm, ...)
center (list): center point at which the regridding is centered to (optional). One of the grid-points is exactely at the given center
refAngle (float): nominal reflection angle of the mirror (float, optional)
Name (string): identifier string of the object (optional)
Raster (float): accuracy of the rastering (step size between points) -- 1 gives four times the points of the input point cloud ==> half spacing
"""
self.printOut("Creating a surface element from a point cloud", self.log.FLOW_TYPE)
if not _GridData :
self.printOut("SciPy-Lib is not available, can't create mirrors from point clouds!", self.log.WARNING_TYPE, False)
return False
self.printOut("Loading a point cloud :", self.log.DEBUG_TYPE, False)
self.printOut("file name : '%s'" % fileName, self.log.DEBUG_TYPE, False)
if len(center) < 3 :
c = np.asarray([0., 0., 0.])
else :
c = np.asarray(center)
# load data file
v = []
z = []
try:
with open(fileName, 'r') as f : # open file for reading
pointNo = 0 # point counter
for line in f : # loop all lines in the file
if ( len(line) == 0 ) or (line[0] == '#') : continue
line = line.strip() # split rows
columns = line.split()
if len(columns) < 3 : continue
if pointNo == 0 : # take first point as center if no center is given
if len(center) < 3 :
c = np.asarray([float(columns[0])*unit, float(columns[1])*unit, float(columns[2])*unit])
else:
c = np.asarray([center[0], center[1], center[1]])
else:
v.append( np.array([ float(columns[0])*unit, float(columns[1])*unit ]) )
z.append( float(columns[2])*unit )
pointNo += 1
except:
self.printOut("File not found!", self.log.ERROR_TYPE, False)
return False
v = np.asarray(v)
z = np.asarray(z)
self.printOut("data file loaded", self.log.DEBUG_TYPE, False)
# rotate mirror to minimize z-range (to keep dz small)
rotStep = 1.
th = rotStep*np.pi/180.
ry = np.asarray([
[ np.cos(th), 0., np.sin(th)],
[ 0., 1., 0.],
[-np.sin(th), 0., np.cos(th)]
])
rx = np.asarray([
[1., 0., 0.],
[0., np.cos(th), -np.sin(th)],
[0., np.sin(th), np.cos(th)]
])
ryInv = np.linalg.inv(ry)
rxInv = np.linalg.inv(rx)
self.printOut("testing rotation around y-axis", self.log.DEBUG_TYPE, False)
# rotate around y-axis
roty = 0.
zSpan = max(z)-min(z)
minzSpan = zSpan + 0.01
while zSpan < minzSpan :
minzSpan = zSpan
for idx, zv in enumerate(z) :
nv = np.dot(ry, [v[idx,0], v[idx,1], zv])
z[idx] = nv[2]
v[idx,0] = nv[0]
v[idx,1] = nv[1]
c = np.dot(ry, c)
roty += rotStep
zSpan = max(z)-min(z)
# de-rotate last step
for idx, zv in enumerate(z) :
nv = np.dot(ryInv, [v[idx,0], v[idx,1], zv])
z[idx] = nv[2]
v[idx,0] = nv[0]
v[idx,1] = nv[1]
c = np.dot(ryInv, c)
roty -= rotStep
self.printOut("found %.1f deg rotation to minimize dz" % roty, self.log.DEBUG_TYPE, False)
self.printOut("testing rotation around x-axis", self.log.DEBUG_TYPE, False)
# rotate around x-axis
rotx = 0
zSpan = max(z)-min(z)
minzSpan = zSpan + 0.01
while zSpan < minzSpan :
minzSpan = zSpan
for idx, zv in enumerate(z) :
nv = np.dot(rx, [v[idx,0], v[idx,1], zv])
z[idx] = nv[2]
v[idx,0] = nv[0]
v[idx,1] = nv[1]
c = np.dot(rx, c)
rotx += rotStep
zSpan = max(z)-min(z)
# de-rotate last step
for idx, zv in enumerate(z) :
nv = np.dot(rxInv, [v[idx,0], v[idx,1], zv])
z[idx] = nv[2]
v[idx,0] = nv[0]
v[idx,1] = nv[1]
c = np.dot(rxInv, c)
rotx -= rotStep
self.printOut("found %.1f deg rotation to minimize dz" % rotx, self.log.DEBUG_TYPE, False)
self.printOut("check out if setup after rotation is convex or concave", self.log.DEBUG_TYPE, False)
# check opening
zout = z[np.argmin(v[:, 0])]
zout += z[np.argmax(v[:, 0])]
zout += z[np.argmin(v[:, 1])]
zout += z[np.argmax(v[:, 1])]
zout = zout / 4.
dmax = np.absolute(max(z) - zout)
dmin = np.absolute(min(z) - zout)
if dmax > dmin : # opening downwards? If so turn around inner and outer behaviour
if inner : inner = False
else : inner = True
self.printOut("found a convex mirror, inverting bool flag 'inner'", self.log.DEBUG_TYPE, False)
else :
self.printOut("found a concave mirror", self.log.DEBUG_TYPE, False)
self.printOut("re-gridding data", self.log.DEBUG_TYPE, False)
# re-grid data
side = int( 2. * np.sqrt(len(v)) / Raster )
minX = min(v[:, 0])
minY = min(v[:, 1])
stepX = (max(v[:, 0]) - min(v[:, 0]) ) / (side-1)
stepY = (max(v[:, 1]) - min(v[:, 1]) ) / (side-1)
sxminus = int((c[0]-min(v[:, 0]))/stepX)
sxplus = int((max(v[:, 0]) - c[0])/stepX)
sideX = sxminus + sxplus + 1
minX = c[0] - sxminus*stepX
maxX = c[0] + sxplus *stepX
syminus = int((c[1]-min(v[:, 1]))/stepY)
syplus = int((max(v[:, 1]) - c[1])/stepY)
sideY = syminus + syplus + 1
minY = c[1] - syminus*stepY
maxY = c[1] + syplus *stepY
# re-grid data
grid_x, grid_y = np.mgrid[minX:maxX:(1j*sideX), minY:maxY:(1j*sideY)]
grid_z = griddata(v, z, (grid_x, grid_y), method='linear')
grid_z.shape = (sideX, sideY)
self.printOut("calculating normal vectors and area element size", self.log.DEBUG_TYPE, False)
# list on non-nan interpolated z-values
inList = np.isnan(grid_z).flatten()
# prepare data array
self.Data = []
self.DataID = (fileName.split("/")[-1]).split(".")[0] # extract filename without extension as
# calculate normal vectors
dist = maxX**2
nC = [0., 0., 0.]
idx = 0
for x in range(0, sideX) : # loop all data points in x
for y in range (0, sideY) : # loop all data points in y
if not inList[idx] : # check if interpolation was successfull
if x == 0 or np.isnan(grid_z[x-1, y]) : # check for grid and mirror edges in x
if not np.isnan(grid_z[x, y]) and not np.isnan(grid_z[x+1, y]) :
dzx = ( grid_z[x+1, y] - grid_z[x, y] )/ stepX # calculate n on starting edge
else : continue
elif x == sideX-1 or np.isnan(grid_z[x+1, y]):
if not np.isnan(grid_z[x, y]) and not np.isnan(grid_z[x-1, y]) :
dzx = ( grid_z[x, y] - grid_z[x-1, y] )/ stepX # calculate n on ending edge
else : continue
else :
if not np.isnan(grid_z[x-1, y]) and not np.isnan(grid_z[x+1, y]) :
dzx = 0.5*( grid_z[x+1, y] - grid_z[x-1, y])/ stepX# calculate n on mirror surface
else : continue
if y == 0 or np.isnan(grid_z[x, y-1]): # check for grid and mirror edges in y
if not np.isnan(grid_z[x, y]) and not np.isnan(grid_z[x, y+1]) :
dzy = ( grid_z[x, y+1] - grid_z[x, y] )/ stepY # calculate n on starting edge
else : continue
elif y == sideY-1 or np.isnan(grid_z[x, y+1]):
if not np.isnan(grid_z[x, y]) and not np.isnan(grid_z[x, y-1]) :
dzy = ( grid_z[x, y] - grid_z[x, y-1] )/ stepY # calculate n on ending edge
else : continue
else :
if not np.isnan(grid_z[x, y-1]) and not np.isnan(grid_z[x, y+1]) :
dzy = 0.5*( grid_z[x, y+1] - grid_z[x, y-1])/ stepY# calculate n on mirror surface
else : continue
# define UFO data point
d = UFODataPoint()
d.v = _UFOVector(minX + x*stepX, minY + y*stepY, grid_z[x,y])
nNorm = np.sqrt(dzx**2+ dzy**2 + 1.)
d.df = stepX*stepY*nNorm
# react on boolean flag 'inner' by inverting the normal vector
if not inner :
d.n = _UFOVector(dzx/nNorm, dzy/nNorm, -1./nNorm)
else :
d.n = _UFOVector(-dzx/nNorm, -dzy/nNorm, 1./nNorm)
self.Data.append(d)
diff = (c[0]-d.v.x)**2 + (c[1]-d.v.y)**2
if diff < dist :
nC = d
dist = diff
idx += 1
self.points = len(self.Data)
self.Header = UFOHeaderPoint()
self.Header.v = nC.v
self.Header.n = nC.n
self.Header.o = _UFOVector(self.Header.n.z, self.Header.n.y, -self.Header.n.x)
self.Data = (UFODataPoint*self.points)(*self.Data)
self.SurfaceType = "PointCloud"
self.FocalLength = None
self.RefAngle = refAngle
self.Name = Name
self.printOut ("got %d mirror data points" % self.points, self.log.DEBUG_TYPE, False)
# back-rotation to original data-point orientation
self.printOut ("rotating back to original position", self.log.DEBUG_TYPE, False)
self.rotX(-rotx)
self.rotY(-roty)
return True
[docs]
def makeLens ( self, f, radius = 10, step = 1., er = 2.25, tanD = 0., Name = "" ) :
"""
Experimental calculation of a thick lens according to Lensmaker's equation.
Parameters:
f (float): focal length [mm]
radius (float): size (radius) of the lens [mm]
step (float): step-size for the discretization [mm]
er (float): permittivity of the lens material
tanD (float): Loss-factor (or tangens delta)
Name (string): identifier string of the object (optional)
Returns:
r2_ufo (UFOData): element with the second lens surface. Simulate first towards the main element and then using the permittivity of the lens towards this second element
"""
self.printOut("Creating a thick Lens from two spherical surfaces", self.log.FLOW_TYPE)
self.printOut("Focal-length : %.2fmm" % f, self.log.DATA_TYPE, False)
self.printOut("Permittivity : %.2f" % er, self.log.DATA_TYPE, False)
self.printOut("tan-Delta : %.2f" % tanD, self.log.DATA_TYPE, False)
self.printOut("Radius : %.2fmm" % radius, self.log.DATA_TYPE, False)
self.printOut("Step-size : %.2fmm" % step, self.log.DATA_TYPE, False)
self.Data = [] # clear data
self.Header = UFOHeaderPoint() # and set data-header
self.Header.v = _UFOVector(0., 0., 0.)
self.Header.n = _UFOVector(0., 0., 1.)
self.Header.o = _UFOVector(1., 0., 0.)
self.RefAngle = 180.
self.Type = _SURFACE # set type to "surface"
offset = 1. # minimum offset between the two surfaces (must be > 0 to avoid calculation errors
n = np.sqrt(er) # calculate diffraction index
if n == 1. : n = 1.01 # and make sure that n-1 is > 0
startAngle = 10.
while True :
R1 = radius/np.sin(np.pi/180*startAngle) # starting radius of the lens (no angle > 10 deg
R2 = 1./(1./R1 - 1./(f*n-f)) # calculate second radius for thin-lens approximation
if np.absolute(R2) < radius*1.2 :
startAngle += 1.
if startAngle > 89. : break
else : break
if startAngle > 89. :return None
p1 = np.cos(np.pi/180.*startAngle)*R1 # calculate lens-thickness
p2 = np.cos((np.arcsin(radius/R2)))*R2
d = R1-p1 + R2-p2 + offset
D = (n-1.)*(1./R1-1./R2+(n-1.)*d/(n*R1*R2))
for stepIter in (10, 1, 0.1, 0.01) : # iteratively calculate R2 of the thick lens by reducing the step-size
while D < 1/f :
R2 += stepIter
p2 = np.cos((np.arcsin(radius/R2)))*R2
d = R1-p1 + R2-p2 + offset
D = (n-1)*(1/R1-1/R2+(n-1)*d/(n*R1*R2))
R2 += stepIter*R2/np.absolute(R2)
p2 = np.cos((np.arcsin(radius/R2)))*R2
d = R1-p1 + R2-p2 + offset
D = (n-1)*(1/R1-1/R2+(n-1)*d/(n*R1*R2))
# create surface R1 as sphere (using the opening angle and the step-size transfered to angles)
r1_angle = np.absolute( np.arcsin(radius/R1)*180./np.pi)
r1_step = np.absolute(np.arcsin(step/p1)*180./np.pi)
self.makeSphere(np.absolute(R1), r1_angle, stepFunction = lambda x : x+r1_step, Name = Name)
if R1 > 0. :
self.rotY(180.)
self.move(0., 0., p1 - R1/np.absolute(R1)*offset)
self.Header.v.z = 0.
self.er = er
# create surface R2 as sphere (using the opening angle and the step-size transfered to angles)
r2_ufo = UFOData(log = self.log)
r2_angle = np.absolute(np.arcsin(radius/R2)*180./np.pi)
r2_step = np.absolute(np.arcsin(step/R2)*180./np.pi)
r2_ufo.makeSphere(np.absolute(R2), r2_angle, stepFunction = lambda x : x+r2_step, Name = Name + " (2)")
if R2 > 0. :
r2_ufo.rotY(180.)
r2_ufo.move(0., 0., p2 - R2/np.absolute(R2)*offset )
r2_ufo.Header.v.z = 0.
self.SurfaceType = "Lens"
self.FocalLength = f
r2_ufo.SurfaceType = "Lens"
r2_ufo.RefAngle = 180.
r2_ufo.FocalLength = f
r2_ufo.er = 1.
r2_ufo.tanD = tanD
return r2_ufo
[docs]
def makeAsphLens ( self, f, radius = 10, step = 1., er = 2.25, tanD = 0., EdgeThickness = 2., Fresnel = None, minThick = 2., Name = "" ) :
"""
Experimental calculation of an aspherical lens.
Parameters:
f (float): focal length [mm]
radius (float): size (radius) of the lens [mm]
step (float): step-size for the discretization [mm]
er (float): permittivity of the lens material
tanD (float): Loss-factor (or tangens delta)
EdgeThickness (float): Thickness at the lens edge [mm]
Fresnel (float or None): Frequency [GHz] at which Fresnel lens steps should be implemented (None == no steps)
Name (string): identifier string of the object (optional)
Returns:
rb_ufo (UFOData): element with the second lens surface. Simulate first towards the main element and then using the permittivity of the lens towards this second element
"""
if isinstance(Fresnel, float) :
LAMBDA = np.sqrt(er)/ (np.sqrt(er)-1)
fresnel = _C/Fresnel *1e-6 *LAMBDA
else :
fresnel = np.inf
self.printOut("Creating an asperical lens", self.log.FLOW_TYPE)
self.printOut("Focal-length : %.2fmm" % f, self.log.DATA_TYPE, False)
self.printOut("Permittivity : %.2f" % er, self.log.DATA_TYPE, False)
self.printOut("tan-Delta : %.2f" % tanD, self.log.DATA_TYPE, False)
self.printOut("Radius : %.2fmm" % radius, self.log.DATA_TYPE, False)
self.printOut("Step-size : %.2fmm" % step, self.log.DATA_TYPE, False)
start = int(radius/step)*step # make sure central point is a data point
self.Data = [] # clear data
self.Header = UFOHeaderPoint() # and set data-header
self.Header.v = _UFOVector( 0, 0, 0)
self.Header.n = _UFOVector( 0, 0, 1)
self.Header.o = _UFOVector(-1, 0, 0)
self.Type = _SURFACE # set type to "surface"
n = np.sqrt(er) # calculate diffraction index
if n == 1. : n = 1.01 # and make sure that n-1 is > 0
# calculate aspherical lense shape iteratively from center (see Wikipedia) along the radius to outside
# reference plane for the lens is the highest point
R = f/(n/(n-1.)-1.) # R0
d = 0. # delta0
z = 0. # z0
contour = [] # new data set
cstep = _STEPLENS/f # step size for the iterative curvature is according to the focal length
for H in np.arange(cstep, radius+3.*cstep, cstep): # loop all data points along radius (central point is already defined)
contour.append(z) # append actual data point
z = d + R - np.sqrt(R*R - H*H) # calculate next data point
gama = np.arctan2(H,f+z) # calculate angle of diffraction
d = n*H/np.sin(gama) # calculate new delta value (first step)
R = np.sqrt( (d - f - z)**2. + H**2.) # calculate new radius of curvature
d = d-R-f # finalize calculation of delta
contour = np.asarray(contour) # transform data array to numpy array
height = np.amax(contour)+EdgeThickness # calculate lens minimum thickness
contour = -contour # and transform lens data to lens curvature
# calculate deviation from curvature data for the normal vector
c1 = contour[1]
c2 = contour[0]
dc = []
for c in contour[1:] :
dc.append(0.5*(c-c1)/cstep)
c1 = c2
c2 = c
dc = np.asarray(dc)
xR = np.arange(0., radius+1.*cstep, cstep)
height = 0.
# calculate lens data points in 3D
for x in np.arange(-start, start, step) : # Loop x-coodinates
for y in np.arange(-start, start, step) : # Loop y coordinates
d = np.sqrt(x*x+y*y) # calculate distance to center
if d > radius : continue # if smaller than lens radius, use data point
alpha = np.arctan2(y, x) # calculate angle with x-axis
dp = UFODataPoint()
# interpolate new data point from curvature data
actz = np.interp(d, xR, contour[:-1])
while actz < -fresnel - minThick: actz += fresnel
if actz < height : height = actz
dp.v = _UFOVector(x, y, actz)
n = np.interp(d, xR, dc)
# and calculate normal vector
n = -1*np.asarray([-n*np.cos(alpha), -n*np.sin(alpha), 1.])
nNorm = np.sqrt(np.dot(n, n))
# and area element size
dp.df = step*step*nNorm
dp.n = _UFOVector(n[0]/nNorm, n[1]/nNorm, n[2]/nNorm)
self.Data.append(dp)
self.points = len(self.Data)
self.Data = (UFODataPoint*self.points)(*self.Data)
self.SurfaceType = "Lens"
self.FocalLength = f
self.RefAngle = 0.
self.Name = Name
rb_ufo = deepcopy(self)
rb_ufo.SurfaceType = "Lens"
for idx in range(0, self.points) :
self.Data[idx].v.z = height-EdgeThickness
self.Data[idx].n = _UFOVector( 0., 0., 1.)
self.Data[idx].df = step*step
self.Header.n = _UFOVector(0., 0., -1.)
self.Header.o = _UFOVector(1., 0., 0.)
self.er = er
rb_ufo.SurfaceType = "Lens"
rb_ufo.FocalLength = f
rb_ufo.Name += " (2)"
rb_ufo.er = 1.
rb_ufo.tanD = tanD
return rb_ufo
#####################################
# plotting
#
# data extraction methods
[docs]
def getCutData ( self, pol = [1., 1., 1.], direction = [0.,90.], axScale = 1., maxDiff = .1, useMax = False) :
"""
Extracts the far-field intensity response of the actual data set along phi direction 'direction'.
Parameters:
pol (list): np.array of floats denoting the polarization vector (will be normalized) used to calculate the intensity
direction (list): of phi-angles with cut-directions or float for a single direction
axScale (float): scaling factor for the x- and y- coordinates (e.g. to have it in opening angles)
maxDiff (float): difference in point distance for using point in cut. (in units of the step-size)
useMax (boolean): indicating if cut passes header coordinates (False) or the field-maximum point (True)
Returns:
tuple (tuple): position (list of np.arrays with the positions), values (list of np.arrays with the cut-data in linear intensity scale)
"""
if isinstance(direction, float) or isinstance(direction, int) :
direction = [direction]
# normalize polarization vector
pol = np.asarray(pol)
pol = pol/ np.sqrt(np.dot(pol, pol))
if self.Type == _SURFACE :
self.printOut("Surface data have no intensity information, exiting!", self.log.ERROR_TYPE)
return None, None
self.printOut("Extracting far field intensity distribution along tracks of different phi-angles", self.log.DEBUG_TYPE)
x = np.zeros(self.points, dtype = float)
y = np.zeros(self.points, dtype = float)
I = np.zeros(self.points, dtype = float)
# extract data
Imax = 0.
idxMax = -1
for idx in range( 0, self.points ) :
x[idx] = self.Data[idx].v.x - self.Header.v.x
y[idx] = self.Data[idx].v.y - self.Header.v.y
E_ = np.array( [ self.Data[idx].Er.x+1j*self.Data[idx].Ei.x, self.Data[idx].Er.y+1j*self.Data[idx].Ei.y, self.Data[idx].Er.z+1j*self.Data[idx].Ei.z ] )
e = np.multiply( E_, pol )
I[idx] = np.absolute( ( 0.5/_ETA )*np.dot( e, e.conj() )*self.Data[idx].df )
Iall = np.absolute( ( 0.5/_ETA )*np.dot( E_, E_.conj() )*self.Data[idx].df )
if Iall > Imax:
Imax = Iall
idxMax = idx
if useMax :
self.printOut("using center position (%.2f, %.2f) at maximum intensity of %.2f" %(x[idxMax], y[idxMax], Imax), self.log.DEBUG_TYPE, False)
x -= x[idxMax]
y -= y[idxMax]
# scale and normalize data
x = x * axScale
y = y * axScale
I = I/Imax
# get unit-cell side size from number of points
unit = 2 * np.amax(np.sqrt(x**2. + y**2))/np.sqrt(len(x))/np.sqrt(2)
# stepList = np.mgrid[ -diff[idxList[-1]]: diff[idxList[-1]]:1j*int( np.sqrt( len(x) ) * 10 ) ]
values = []
position = []
for idx, d in enumerate(direction) :
vals = []
posis = []
cos = np.cos( np.pi*d/180. )
sin = np.sin( np.pi*d/180. )
bIdx = np.asarray(np.absolute(x*sin+y*cos) < unit*maxDiff) # rotate y-coordinates and take bool list of those which near zero
position.append( deepcopy((x*cos-y*sin)[bIdx]) ) # rotate x-coordinates and take the ones from bool-list
values.append ( deepcopy(I[bIdx]) ) # take belonging intensities via bool-list
return position, values
[docs]
def getCutDataInter ( self, pol = [1., 1., 1.], direction = [0.,90.], axScale = 1., useMax = False) :
"""
Extracts the far-field intensity response of the actual data set along phi direction 'direction'.
Parameters :
pol (list): np.array of floats denoting the polarization vector (will be normalized) used to calculate the intensity
direction (list): of phi-angles with cut-directions or float for a single direction
axScale (float): scaling factor for the x- and y- coordinates (e.g. to have it in opening angles)
useMax (boolean): indicating if cut passes header coordinates (False) or the field-maximum point (True)
Returns:
tuple (tuple): position (list of np.arrays with the positions), values (list of np.arrays with the cut-data in linear intensity scale)
"""
if isinstance(direction, float) or isinstance(direction, int) :
direction = [direction]
# normalize polarization vector
pol = np.asarray(pol)
pol = pol/ np.sqrt(np.dot(pol, pol))
if self.Type == _SURFACE :
self.printOut("Surface data have no intensity information, exiting!", self.log.ERROR_TYPE)
return None, None
self.printOut("Extracting far field intensity distribution along tracks of different phi-angles", self.log.DEBUG_TYPE)
x = np.zeros(self.points, dtype = float)
y = np.zeros(self.points, dtype = float)
I = np.zeros(self.points, dtype = float)
Imax = 0.
idxMax = -1
for idx in range( 0, self.points ) :
x[idx] = self.Data[idx].v.x - self.Header.v.x
y[idx] = self.Data[idx].v.y - self.Header.v.y
E_ = np.array( [ self.Data[idx].Er.x+1j*self.Data[idx].Ei.x, self.Data[idx].Er.y+1j*self.Data[idx].Ei.y, self.Data[idx].Er.z+1j*self.Data[idx].Ei.z ] )
e = np.multiply( E_, pol )
I[idx] = np.absolute( ( 0.5/_ETA )*np.dot( e, e.conj() )*self.Data[idx].df )
Iall = np.absolute( ( 0.5/_ETA )*np.dot( E_, E_.conj() )*self.Data[idx].df )
if Iall > Imax:
Imax = Iall
idxMax = idx
if useMax :
x -= ( self.Data[idxMax].v.x - self.Header.v.x )
y -= ( self.Data[idxMax].v.y - self.Header.v.y )
x = x * axScale
y = y * axScale
I = I/Imax
# get an approximation of the median point-difference
diff = np.sqrt(x**2. + y**2)
idxList = np.argsort(diff)
unit = 2.*diff[idxList[-1]]/np.sqrt(len(x))
stepList = np.mgrid[ -diff[idxList[-1]]: diff[idxList[-1]]:1j*int( np.sqrt( len(x)*2 ) ) ]
values = []
position = []
for idx, d in enumerate(direction) :
vals = []
posis = []
for pos in stepList :
xpos = np.cos( np.pi*d/180. )*pos
ypos = np.sin( np.pi*d/180. )*pos
if (xpos >= np.min(x)) and (xpos <= np.max(x)) and (ypos >= np.min(y)) and (ypos <= np.max(y)) :
# calculate differences to available points
diff = (x-xpos)**2. + (y-ypos)**2
idxList = np.argsort(diff)
# check if it is close to one point
if (diff[idxList[0]] < unit/100000.) :
value = I[idxList[0]] # then use this point
else :
value = 0. # otherwise interpolate between next 8 neighbors
norm = 0.
for i in range(0, 4):
value += I[idxList[i]]/diff[idxList[i]]
norm += 1./diff[idxList[i]]
value = value/norm
vals.append( value )
posis.append( pos )
position.append( np.asarray( posis ) )
values.append ( np.asarray( vals ) )
return position, values
# plotting methods
[docs]
def plotIntensity ( self, figname = "", show = False, pol = [1., 1., 1.], dB = False, Field = False, returnFIG = False, axScale = 1., axUnit = "mm") :
"""
Plot the far-field intensity response of the actual data set.
Parameters:
figname (string) : Filename 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)
pol (list) : polarization vector (will be normalized) used to calculate the intensity
dB (boolean) : indicating if scale is in dB pr linear
Field (boolean) : indicating if the intensity or the field-amplitude is plotted
returnFIGt (boolean) : indicating if figure handle will be returned after plotting (figname and show have no effect) or not
axScalet (float) : scaling factor for the x- and y- coordinates (e.g. to have it in opening angles)
axInit (string) : unit string for the x- and y- axis
Returns:
tuple (tuple): (fig, ax) figure description and axis (if returnFIG is set) or None, None
"""
self.printOut("Plotting far field intensity distribution", self.log.FLOW_TYPE)
if not _GridData :
self.printOut("SciPy-Lib is not available, can't create the plot!", self.log.WARNING_TYPE, False)
return None, None
# normalize polarization vector
pol = np.asarray(pol)
pol = pol/ np.sqrt(np.dot(pol, pol))
if self.Type == _SURFACE :
self.printOut("Surface data have no intensity information, exiting!", self.log.ERROR_TYPE, False)
return
x = np.zeros(self.points, dtype = float)
y = np.zeros(self.points, dtype = float)
I = np.zeros(self.points, dtype = float)
Imax = 0.
for idx in range(0, self.points) :
x[idx] = self.Data[idx].v.x
y[idx] = self.Data[idx].v.y
E_ = np.array( [self.Data[idx].Er.x + 1j * self.Data[idx].Ei.x, self.Data[idx].Er.y + 1j * self.Data[idx].Ei.y, self.Data[idx].Er.z + 1j * self.Data[idx].Ei.z] )
e = np.multiply( E_, pol )
I[idx] = np.absolute( ( 0.5 / _ETA ) * np.dot( e, e.conj() ) * self.Data[idx].df )
Iall = np.absolute( ( 0.5 / _ETA ) * np.dot( E_, E_.conj() ) * self.Data[idx].df )
if Iall > Imax: Imax = Iall
x = x * axScale
y = y * axScale
I = I/Imax
if np.amax(np.absolute(I)) < 0.5 :
levelsdB = [-20, -10, -3]
levelsL = [0.001, 0.01, 0.1]
else :
levelsdB = [-13, -3, -1]
levelsL = [0.1, 0.5, 0.95]
label = "intensity"
if Field:
I = np.sqrt(I)
label = "field amplitude"
if dB :
I = np.log10(I)*10.
levels = levelsdB
else :
levels = levelsL
# re-grid data
xi, yi = np.mgrid[min(x):max(x):100j, min(y):max(y):100j]
zi = griddata((x, y), I, (xi, yi), method = 'linear')
fig, ax = plt.subplots()
cf = ax.pcolormesh(xi, yi, zi, cmap = plt.get_cmap('rainbow'))
if dB :
ax.contour(xi, yi, zi, levels=levels, linewidths = 0.5, colors = 'k')
ax.set_title ("normalized %s [dB],\n polarization vector: (%.1f, %.1f, %.1f),\n(contours are at %.1f, %.1f, and %.1f)"
% (label, pol[0], pol[1], pol[2], levels[0], levels[1], levels[2]))
else :
ax.contour(xi, yi, zi, levels=levels, linewidths = 0.5, colors = 'k')
ax.set_title ("normalized %s,\n polarization selection vector: (%.1f, %.1f, %.1f),\n(contours are at %.2f, %.2f, and %.2f)"
% (label, pol[0], pol[1], pol[2], levels[0], levels[1], levels[2]))
ax.set_xlabel("position [%s]" % axUnit)
ax.set_ylabel("position [%s]" % axUnit)
fig.colorbar(cf)
ax.set_aspect('equal')
if returnFIG : return fig, ax
if figname != "" :
fig.savefig(figname)
if show :
plt.show()
else :
plt.show()
plt.close()
return None, None
[docs]
def plotCut ( self, figname = "", show = False, pol = [1., 1., 1.], direction = [0.,90.], dB = False, Field = False, axScale = 1., axUnit = "mm") :
"""
Plot the far-field intensity response of the actual data set along phi direction direction.
Parameters:
figname (string): Filename 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)
pol (list): polarization vector (will be normalized) used to calculate the intensity
direction (list): of phi-angles with cut-directions or float for a single direction
dB (boolean): indicating if scale is in dB pr linear
Field (boolean): indicating if the intensity or the field-amplitude is plotted
axScalet (float): scaling factor for the x- and y- coordinates (e.g. to have it in opening angles)
axUnit (string): unit string for the x- and y- axis
Returns:
tuple (tuple): (fig, ax) figure description and axis (if returnFIG is set) or None, None
"""
self.printOut("Plotting far field intensity distribution along phi-angles", self.log.FLOW_TYPE)
if self.Type == _SURFACE :
self.printOut("Surface data have no intensity information, exiting!", self.log.ERROR_TYPE, False)
return None, None
# extract data
position, values = self.getCutData ( pol, direction, axScale )
# get line-colors
colors = getColors(len(direction))
# initialize plot
label = "intensity"
if Field:
for idx in range(0, len(values)) :
values[idx] = np.sqrt(values[idx])
label = "amplitude"
if dB :
for idx in range(0, len(values)) :
values[idx] = np.log10(values[idx])*10.
fig, ax = plt.subplots()
if dB :
ax.set_title ("normalized field %s [dB],\n polarization vector: (%.1f, %.1f, %.1f)"
% ( label, pol[0], pol[1], pol[2] ) )
ax.set_ylabel("%s [dB] " % label )
else :
ax.set_title ("normalized field %s,\n polarization selection vector: (%.1f, %.1f, %.1f)"
% ( label, pol[0], pol[1], pol[2] ) )
ax.set_ylabel("%s [au] " % label )
ax.set_xlabel("position [%s]" % axUnit)
# plot all directions
for idx in range(0, len(values)) :
ax.plot(position[idx], values[idx], color = colors[idx], linestyle = 'solid', label = "angle = %.1f degree" % direction[idx])
# plot legend
fig.legend(loc='upper left', bbox_to_anchor=(0.65, 0.85))
# check for saving and showing the image
if figname != "" :
fig.savefig(figname)
if show :
plt.show()
else :
plt.show()
plt.close()
return None, None
[docs]
def plotPhase ( self, figname = "", show = False, pol = [1., 1., 1.], returnFIG = False) :
"""
Plot the far-field phase response of the actual data set.
Parameters:
figname (string): Filename 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)
pol (list): polarization vector (will be normalized) used to calculate the intensity
returnFIGt (boolean): indicating if figure handle will be returned after plotting (figname and show have no effect) or not
Returns:
tuple (tuple): (fig, ax) figure description and axis (if returnFIG is set) or None, None
"""
self.printOut("Plotting far field intensity distribution", self.log.FLOW_TYPE)
if not _GridData :
self.printOut("SciPy-Lib is not available, can't create the plot!", self.log.WARNING_TYPE, False)
return None, None
# normalize polarization vector
pol = np.asarray(pol)
pol = pol/ np.sqrt(np.dot(pol, pol))
if self.Type == _SURFACE :
self.printOut("Surface data have no intensity information, exiting!", self.log.ERROR_TYPE, False)
return
x = np.zeros(self.points, dtype = float)
y = np.zeros(self.points, dtype = float)
A = np.zeros(self.points, dtype = float)
c = 0.
d = 1e10
for idx in range(0, self.points) :
x[idx] = self.Data[idx].v.x
y[idx] = self.Data[idx].v.y
E_ = np.array([self.Data[idx].Er.x+1j*self.Data[idx].Ei.x, self.Data[idx].Er.y+1j*self.Data[idx].Ei.y, self.Data[idx].Er.z+1j*self.Data[idx].Ei.z])
e = np.dot(E_, pol)
A[idx] = np.angle(e, deg = True)
dact = self.Data[idx].v.x**2 + self.Data[idx].v.y**2
if dact < d :
d = dact
c = np.angle(e, deg = True)
A = A - c
for ind in range(0, len(A)) :
while A[ind] > 180. : A[ind] -= 180.
while A[ind] < -180. : A[ind] += 180.
# re-grid data
xi, yi = np.mgrid[min(x):max(x):100j, min(y):max(y):100j]
zi = griddata((x, y), A, (xi, yi), method='linear')
fig, ax = plt.subplots()
cf = ax.pcolormesh(xi, yi, zi, cmap = plt.get_cmap('rainbow'))
ax.contour(xi, yi, zi, levels=[-180, -90, -10, 0, 10, 90, 180], linewidths = 0.5, colors = 'k')
ax.set_title ("Phase of the field [deg],\n polarization selection vector: (%.1f, %.1f, %.1f)),\n countours (-180, -90, -10, 0, 10, -90, 180)"
% (pol[0], pol[1], pol[2]))
ax.set_xlabel("position [mm]")
ax.set_ylabel("position [mm]")
fig.colorbar(cf)
ax.set_aspect('equal')
if returnFIG : return fig, ax
if figname != "" :
fig.savefig(figname)
if show :
plt.show()
else :
plt.show()
plt.close()
return None, None
[docs]
def plotMirror ( self, figname = "", show = False) :
"""
Plot the three 2d projections of the actual mirror data.
Parameters:
figname (string): Filename 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 mirror structure", self.log.FLOW_TYPE)
x = np.zeros(self.points, dtype = float)
y = np.zeros(self.points, dtype = float)
z = np.zeros(self.points, dtype = float)
for idx in range(0, self.points) :
x[idx] = self.Data[idx].v.x
y[idx] = self.Data[idx].v.y
z[idx] = self.Data[idx].v.z
fig, ax = plt.subplots(3)
ax[0].plot(x, y, 'r.')
ax[0].set_aspect('equal')
ax[0].set_title ("x-y projection")
ax[1].plot(x, z, 'r.')
ax[1].set_aspect('equal')
ax[1].set_title ("x-z projection")
ax[2].plot(y, z, 'r.')
ax[2].set_aspect('equal')
ax[2].set_title ("y-z projection")
ax[2].set_xlabel("position [mm]")
ax[2].set_ylabel("position [mm]")
fig.suptitle("Projections of the 3d-mirror structure")
fig.tight_layout(h_pad = 2. )
plt.subplots_adjust(top = 0.85)
if figname != "" :
fig.savefig(figname)
if show :
plt.show()
else :
plt.show()
plt.close()
#####################################
# coordination transformations
[docs]
def getRot ( self ) :
"""
Calculates the rotation-angles to get the header normal vector along the positive z-axis.
"""
self.printOut("Getting rotation angles", self.log.FLOW_TYPE)
if len(self.Data) < 2 :
self.printOut("no or not enough data available -> exiting", self.log.ERROR_TYPE, False)
return 0, 0, 0
n = np.asarray([self.Header.n.x, self.Header.n.y, self.Header.n.z])
aX = np.arctan2(n[1], n[2])
aY = -np.arctan2(n[0], np.sqrt(n[1]**2+n[2]**2))
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)]
])
# aZ : Header.o must be oriented in x-z plane
dnew = np.dot( ry, np.dot( rx, np.asarray( [self.Header.o.x, self.Header.o.y, self.Header.o.z] ) ) )
aZ = -np.arctan2(dnew[1], dnew[0])
return aX*180./np.pi, aY*180./np.pi, aZ*180./np.pi
[docs]
def resetPosition ( self ) :
"""
Removes offset and all rotations (makes sure normal vector of the header point is pointing to (0,0,1)).
"""
# remove offsets
v = self.getCenter()
self.move(-v[0], -v[1], -v[2])
# rotate back
rx, ry, rz = self.getRot()
self.rotX(rx)
self.rotY(ry)
self.rotZ(rz)
[docs]
def getCenter ( self ) :
"""
Get the nominal center of the part according to its heder information.
"""
return np.asarray([self.Header.v.x, self.Header.v.y, self.Header.v.z])
[docs]
def rotX ( self, angle = 0.) :
"""
Rotate UFO data around the x-axis.
Parameters:
angle (float): rotation angle in degree
"""
self.printOut("Rotating around x-axis by %.2f degree" % angle, self.log.DEBUG_TYPE)
th = angle*np.pi/180.
rx = np.asarray([
[1., 0., 0.],
[0., np.cos(th), -np.sin(th)],
[0., np.sin(th), np.cos(th)]
])
self.transform(rx)
[docs]
def rotY ( self, angle = 0.) :
"""
Rotate UFO data around the y-axis.
Parameters:
angle (float): rotation angle in degree
"""
self.printOut("Rotating around y-axis by %.2f degree" % angle, self.log.DEBUG_TYPE)
th = angle*np.pi/180.
ry = np.asarray([
[ np.cos(th), 0., np.sin(th)],
[ 0., 1., 0.],
[-np.sin(th), 0., np.cos(th)]
])
self.transform(ry)
[docs]
def rotZ ( self, angle = 0.) :
"""
Rotate UFO data around the y-axis.
Parameters:
angle (float): rotation angle in degree
"""
self.printOut("Rotating around z-axis by %.2f degree" % angle, self.log.DEBUG_TYPE)
th = angle*np.pi/180.
rz = np.asarray([
[np.cos(th), -np.sin(th), 0.],
[np.sin(th), np.cos(th), 0.],
[ 0., 0., 1.]
])
self.transform(rz)
[docs]
def move ( self, x = 0., y = 0., z = 0., absolute = False ) :
"""
Move the data by the vextor (x,y,z) (all in mm).
"""
# overload the x-parameter and allow vectors instead of three independent numbers
try :
l = len( x )
except :
l = 0
if ( l == 3 ) and ( y == 0 ) and ( z == 0 ):
y = x[1]
z = x[2]
x = x[0]
l = 0
if l > 0 :
self.printOut("Cannot interprete vector input for moving , exiting!", self.log.ERROR_TYPE)
self.printOut("Moving part by (%.2fmm, %.2fmm, %.2fmm)" % (x, y, z), self.log.DEBUG_TYPE)
if absolute:
self.moveZero()
self.Header.v.x += x
self.Header.v.y += y
self.Header.v.z += z
self.Data, self.points = _cMove(self.Data, self.points, _UFOVector(x,y,z))
[docs]
def moveZero ( self ):
self.move( -self.Header.v.x, -self.Header.v.y, -self.Header.v.z )
#####################################
# fit operations
[docs]
def fitBeam ( self, width = 100. ) :
"""
Fit a 2D Gaussian to the intensity distribution. A simple round gaussian with 1/e width 'width' is used.
Parameters:
width (float) : starting value (expectation value) of the 1/e width of the intensity distribution
Returns:
amplitude (float) W
1/e-width (float) [mm]
x-postion (float) [mm]
y-postion (float) [mm]
"""
self.printOut("Fitting a gaussian beam to the intensity distribution" , self.log.FLOW_TYPE)
self.printOut(" using a 1/e width of %.2f mm as starting value" % ( width ), self.log.DATA_TYPE)
ufo = deepcopy( self )
ufo.resetPosition()
x = []
I = []
maxD = ufo.Data[0]
maxSize = ufo.Data[0]
for d in ufo.Data :
x.append([ d.v.x, d.v.y ])
I.append( d.I )
if d.I > maxD.I:
maxD = d
if d.v.x > maxSize.v.x :
maxSize = d
x = np.asarray( x ).T / maxSize.v.x
I = np.asarray( I ) / maxD.I
def _gaussFunc( xy, *Pin ) :
x, y = xy
A = float(Pin[0])
w = float(Pin[1])
xo = float(Pin[2])
yo = float(Pin[3])
return A * np.exp( -( ( x - xo )**2 + ( y - yo )**2 ) / ( w * w ) )
Pstart = [1., width/maxSize.v.x, int(maxD.v.x/maxSize.v.x), int(maxD.v.y/maxSize.v.x) ]
if not _GridData :
return Pstart
popt, pcov = curve_fit( _gaussFunc, x, I, p0 = Pstart, ftol = 1e-12, method = 'lm' )
popt[0] *= maxD.I
popt[1] *= maxSize.v.x
popt[2] *= maxSize.v.x
popt[3] *= maxSize.v.x
return popt
#####################################
# data operations
def _modField ( self, factor = 1. + 1j, onlyE = True, JRef = True, JTrans = True) :
"""
Modifies the internal field data by multiplying it with the complex factor "factor".
Parameters:
factor (complex float): factor to multipy fields and currents with (intensity is multiplied with the square of the absolute value)
onlyE (boolean): indicating if True, that the currents are not multiplied
JRef (boolean): indicating if the reflected currents are modified
JTrans (boolean): indicating if the transmitted currents are modified
"""
s = "and currents"
if onlyE :
s = ""
self.printOut("Modifying field %s by factor %.2f" % ( s, factor.real ), self.log.DEBUG_TYPE)
for d in self.Data :
if not onlyE :
if JRef :
d.Jr.x *= factor
d.Jr.y *= factor
d.Jr.z *= factor
d.Ji.x *= factor
d.Ji.y *= factor
d.Ji.z *= factor
if JTrans :
d.Jtr.x *= factor
d.Jtr.y *= factor
d.Jtr.z *= factor
d.Jti.x *= factor
d.Jti.y *= factor
d.Jti.z *= factor
d.Er.x *= factor
d.Er.y *= factor
d.Er.z *= factor
d.Ei.x *= factor
d.Ei.y *= factor
d.Ei.z *= factor
d.I *= np.absolute( factor ) ** 2.
def _phaseShift ( self, shift = 180., onlyJ = True, JRef = True, JTrans = True) :
"""
Phase shift of the internal fields
Parameters:
onlyE (boolean): indicating if True, that only the currents undergo a phase shift
JRef (boolean): indicating if the reflected currents are modified
JTrans (boolean): indicating if the transmitted currents are modified
"""
s = "and E-fields"
if onlyJ :
s = "only"
self.printOut("introducing phase shift of %.2f deg for currents %s" % ( shift, s ), self.log.DEBUG_TYPE)
shift = np.pi * shift / 180.
factor = np.exp( 1j * shift )
for d in self.Data :
# reflected currents
if JRef :
jxr = ( d.Jr.x + 1j*d.Ji.x ) * factor
jyr = ( d.Jr.y + 1j*d.Ji.y ) * factor
jzr = ( d.Jr.z + 1j*d.Ji.z ) * factor
d.Jr.x = jxr.real
d.Jr.y = jyr.real
d.Jr.z = jzr.real
d.Ji.x = jxr.imag
d.Ji.y = jyr.imag
d.Ji.z = jzr.imag
# transmitted currents
if JTrans :
jxt = ( d.Jtr.x + 1j*d.Jti.x ) * factor
jyt = ( d.Jtr.y + 1j*d.Jti.y ) * factor
jzt = ( d.Jtr.z + 1j*d.Jti.z ) * factor
d.Jtr.x = jxt.real
d.Jtr.y = jyt.real
d.Jtr.z = jzt.real
d.Jti.x = jxt.imag
d.Jti.y = jyt.imag
d.Jti.z = jzt.imag
# E-field
if not onlyJ :
ex = ( d.Er.x + 1j*d.Ei.x ) * factor
ey = ( d.Er.y + 1j*d.Ei.y ) * factor
ez = ( d.Er.z + 1j*d.Ei.z ) * factor
d.Er.x = ex.real
d.Er.y = ey.real
d.Er.z = ez.real
d.Ei.x = ex.imag
d.Ei.y = ey.imag
d.Ei.z = ez.imag
[docs]
def applyLosses ( self, lossFactor, lossTemperature, phCenter = [], verbose = True ) :
"""
Applies losses to the existing fields. It calculates the loss in respect to the loss-temperature and the isometric temperature from constants.
Parameters:
lossFactor (float): 0 to 1., denoting the losses in the field during the trace.
lossTemperature (float): temperature [K] at which losses are dumped
phaseCenter (list): vector denoting the face center of the reflected signal. must be in reflection direction, should be in distance of the waist and at least 100 times the wavelength away.
verbose (boolean): quiet opertion if True
Returns:
tuple (tuple): ufo (UFODataClass with loss-power), ufoT (UFODataClass with temperature-normalized loss-power)
"""
# implement losses
if verbose:
self.printOut("Applying losses of %.2f %% which are terminated at a temperature of %.2f K :" %(lossFactor*100, lossTemperature), self.log.FLOW_TYPE)
ufo = deepcopy( self )
self._modField ( np.sqrt(1. - lossFactor), onlyE = False, JRef = True, JTrans = True ) # apply losses to all outputs (curents)
ufo._modField ( np.sqrt( lossFactor ), onlyE = False, JRef = True, JTrans = True ) # create loss-object with remaining power (Transmitted and reflected)
if len(phCenter) == 3 : # is a phase center given to calculate a power sphere from?
pR = 0. # if so get reflected and transmitted loss-power
pT = 0.
for d in ufo.Data :
jR = np.asarray( [d.Jr.x + 1j*d.Ji.x, d.Jr.y + 1j*d.Ji.y, d.Jr.z + 1j*d.Ji.z] )
jT = np.asarray( [d.Jtr.x + 1j*d.Jti.x, d.Jtr.y + 1j*d.Jti.y, d.Jtr.z + 1j*d.Jti.z] )
pR += np.dot( jR, jR.conj() )
pT += np.dot( jT, jT.conj() )
if pR != 0 : # check if reflected loss-power is present at all
lossCorretion = ( ( pT + pR ) / pR ).real # calculate factor to get full power from reflected power
ufo._modField ( lossCorretion, onlyE = False ) # and correct reflected fields for the transmitted part of the power
else : # else use transmitted power alone
for d in ufo.Data : # copy conj transmitted currents to real currents and use them in the following
d.Jr.x = d.Jtr.x
d.Jr.y = d.Jtr.y
d.Jr.z = d.Jtr.z
d.Ji.x = -d.Jti.x
d.Ji.y = -d.Jti.y
d.Ji.z = -d.Jti.z
# create a sphere with center at the mirror origin, radius is touching the phase center and normal vector is aligned in direction towards phase center
phCenter = np.asarray( phCenter )
n = np.asarray( [ufo.Header.v.x, ufo.Header.v.y, ufo.Header.v.z] ) - phCenter
r = np.sqrt( np.dot( n, n ) ) # get sphere radius (length of the distance to phase center)
n = n / r # get normal vector (normalized connection vector between mirror and phase center)
r = 1000. * self.Wavelength
if verbose :
self.printOut("using a phase center of absolute (%.2f mm, %.2fmm, %.2fmm)" % (phCenter[0], phCenter[1], phCenter[2]), self.log.DEBUG_TYPE, False )
self.printOut("resulting in a light direction of (%.2f, %.2f, %.2f)" % (n[0], n[1], n[2]), self.log.DEBUG_TYPE, False )
# create UFO sphere (4 pi )
target = UFOData()
# calculate maximum defraction angle
Aeff = self.getAeff( n )
D = 2. * np.sqrt( Aeff * 4. /np.pi ) # calculate effective diameter of the surface
raster = np.sqrt( Aeff / self.points ) # effective raster
maxTheta = 0.8 * np.arctan2( self.Wavelength, raster ) * 180. / np.pi
# define step-function to create sphere
def stepFunc ( x ) :
nStep = 2. # normal step-size
angleDense = 4 * self.Wavelength / D * 180. / np.pi # step size in 4 FWHM area
angleMedi = 10 * self.Wavelength / D * 180. / np.pi # step size in 10 FWHM area
angleLarge = 25 * self.Wavelength / D * 180. / np.pi # step size in 10 FWHM area
stepDense = angleDense / 40.
stepMedi = angleMedi / 50.
stepLarge = angleLarge / 50.
if np.abs( x ) < angleDense : return x + stepDense
elif np.abs( x ) < angleMedi : return x + stepMedi
elif np.abs( x ) < angleLarge : return x + stepLarge
return x + nStep
target.makeSphere( radius = r, thetaMax = maxTheta, stepFunction = stepFunc )
target.Header.n = _UFOVector( -n[0], -n[1], -n[2] ) # get rotation angles via class function, therefore modify header entry and set it back later
aX, aY, aZ = target.getRot()
target.Header.n = _UFOVector(0., 0., -1.) # restore original normal vector
target.rotZ ( -aZ ) # rotate normal vector along connection vector
target.rotY ( -aY )
target.rotX ( -aX )
aX, aY, aZ = ufo.getRot( ) # get rotation angles of self surface
target.rotX ( aX ) # and rotate sphere by these angles, all data points with negative z-value are behind the mirror
target.rotY ( aY )
target.rotZ ( aZ )
# filter points to use only points in front (positive z-value)
data = []
for d in target.Data :
if d.v.z > r / 20. :
data.append(d)
target.points = len( data ) # set new data points to class
target.Data = ( UFODataPoint * target.points )( *data )
target.rotZ ( -aZ ) # rotate normal vector along connection vector
target.rotY ( -aY )
target.rotX ( -aX )
target.move ( phCenter[0], phCenter[1], phCenter[2] ) # move sphere center on mirror-center (and sphere mirror origin to phase center)
target.Name = ufo.Name + "-Power-Sphere"
ufo = ufo.simulate( target, self.Wavelength ) # simulate e-field on the sphere
ufoT = deepcopy( ufo ) # create temperaure weighted loss-class
ufoT._modField ( np.sqrt( lossTemperature / _Tiso ), onlyE = False )
if verbose:
self.printOut("Power in loss class : %.2f" % ( ufo.getPower()[0] ), self.log.DEBUG_TYPE, False)
self.printOut("Power on mirror : %.2f" % ( self.getPower(phCenter = phCenter )[0] ), self.log.DEBUG_TYPE, False)
return ufo, ufoT
[docs]
def getAeff ( self, d = [0., 0., 1.] ) :
"""
Calculates the effective area of the ufo surface along the projection of d (not for intensity files).
Parameters:
d (list): ???
Returns:
Aeff (float): effective area
"""
d = np.array(d)
dn = np.sqrt(d[0]*d[0]+d[1]*d[1]+d[2]*d[2])
d = d / dn
Aeff = 0.
for idx in range(0, self.points) :
if ( self.Type == _INTENSITY ) :
Aeff += self.Data[idx].df
else :
Aeff += self.Data[idx].df * np.absolute( np.dot( d, np.asarray( [self.Data[idx].n.x, self.Data[idx].n.y, self.Data[idx].n.z] ) ) )
self.printOut("Total surface area of Dataset %s is %.2f mm^2" % (self.DataID, Aeff), self.log.DEBUG_TYPE)
return Aeff
[docs]
def correlate ( self, ufo2) :
"""
Calculates overlap integral between this ufo-class and class given in ufo2.
Parameters:
ufo2 (UFODataClass): UFODataClass with calculated fields. The function will calculate the overlap integral. The both classes (self and this) must share a common point cloud.
Returns:
tuple (tuple): two floats, overlap integral between e-fields and normalization factor
"""
self.printOut("Correlating Data :", self.log.FLOW_TYPE)
self.printOut("Data-set 1 : Privat data (Data-ID : '%s')" % self.DataID, self.log.DEBUG_TYPE, False)
self.printOut("Data-set 2 : Data-ID '%s'" % ufo2.DataID, self.log.DEBUG_TYPE, False)
# check if data ranges are comatible
if self.points != ufo2.points :
self.printOut("Data-Sets have different length (%d and %d points), correlation not possible!" % (self.points, ufo2.points), self.log.ERROR_TYPE, False)
return -1., -1.
# check if file type is correct
if ( self.Type == _SURFACE) or (ufo2.Type == _SURFACE) :
self.printOut("Data-Sets have wrong type (at least one is a UFO-Surface)", self.log.ERROR_TYPE, False)
return -1., -1.
# reset all variables
i1 = 0. # overall intensity of self.Data
i2 = 0. # overall intensity of ufo2.Data
ol = 0.+1j*0.
for idx in range(0, self.points) :
d1 = self.Data[idx]
d2 = ufo2.Data[idx]
E1 = np.array([d1.Er.x+1j*d1.Ei.x, d1.Er.y+1j*d1.Ei.y, d1.Er.z+1j*d1.Ei.z])
E2 = np.array([d2.Er.x+1j*d2.Ei.x, d2.Er.y+1j*d2.Ei.y, d2.Er.z+1j*d2.Ei.z])
i1 += (0.5/_ETA) * np.dot(E1, E1.conj()) * d1.df
i2 += (0.5/_ETA) * np.dot(E2, E2.conj()) * d2.df
ol += (0.5/_ETA) * np.dot(E1, E2.conj()) * d1.df
self.printOut("Correlation factor = %.3f %+.3fi" %(ol.real, ol.imag), self.log.DEBUG_TYPE, False)
self.printOut("Normalisazion = %.3f" %(np.sqrt(i1*i2).real), self.log.DEBUG_TYPE, False)
return np.absolute(ol), np.sqrt(i1*i2).real
[docs]
def setI0 ( self, I0 = -1, setOnly = True ) :
"""
Sets the internal excitation current to either the given value (if > 0) or calculates it from the overall power.
Parameters:
I0 (float): new exitation current [A]; if negative curent is calculated from actual available power in the field
setOnly (bool): if True, only current will be set, if False, field strength is adapted to new current
"""
if setOnly and ( I0 > 0 ):
self.I0 = I0
return
power = self.getPower()
if power != None :
self.I0 = np.sqrt(2.*power[0]*1e-6/_Impedance)
if I0 > 0 :
factor = ( I0 / self.I0 )
self._modField ( factor = factor, onlyE = False, JRef = True, JTrans = True)
self.I0 = I0
[docs]
def getPower ( self, radius = -1, polVector = [1., 1., 1.], phCenter = []) :
"""
Calculates the power in the field at the surface.
Attention, this only works if the surface is a sphere!!
A flat surface at waist position is also fine (spherical radius of infinity, flat phase distribution).
Giving the phase center, and therewith a direction of propagation and a distance for calculation, power is calculated correctly also on curved surfaces (timeconsuming).
Parameters:
radius (float): if > 0, use only points within the circle given by radius
polVector (list): list or array with the components of a polarization vector (will be normalized) on which the fields are projected.
phaseCenter (list): vector denoting the phase center of the reflected signal. Must be in reflection direction, should be in distance of the waist and at least 100 times the wavelength away.
"""
self.printOut("calculating the overall power", self.log.FLOW_TYPE)
polVector = np.asarray(polVector)
# check if file type is correct
if ( self.Type == _SURFACE):
self.printOut("Data-Set has wrong type (UFO-Surface)", self.log.ERROR_TYPE, False)
return 0., 0.
if len(phCenter) == 3 :
use = deepcopy( self )
if radius > 0. :
data = []
headerV = use.Header.v
for d in use.Data :
if np.sqrt( (d.v.x - headerV.x)**2 + (d.v.y - headerV.y)**2 + (d.v.z - headerV.z)**2 ) < radius :
data.append(d)
use.points = len( data )
use.Data = ( UFODataPoint * use.points )( *data )
use, useT = use.applyLosses ( 1., 0., phCenter = phCenter, verbose = False )
return use.getPower()
else :
# reset all variables
i = 0. # overall intensity of self.Data
F = 0. # overall area
headerV = self.Header.v
for c in range(0, self.points) :
d = self.Data[c]
r = np.sqrt( (d.v.x - headerV.x)**2 + (d.v.y - headerV.y)**2 + (d.v.z - headerV.z)**2)
if (radius < 0) or ( r < radius) :
E = np.array([d.Er.x+1j*d.Ei.x, d.Er.y+1j*d.Ei.y, d.Er.z+1j*d.Ei.z])
E = np.multiply(E, polVector)
F += d.df
i += (np.sqrt(self.er) * 0.5/_ETA) * np.dot(E , E.conj()) * d.df
self.printOut(" power = %.3f" %(i.real), self.log.DEBUG_TYPE, False)
return i.real, F
[docs]
def getPolarization ( self ) :
self.printOut("Calculating the linear polarization vector", self.log.FLOW_TYPE)
# check if file type is correct
if ( self.Type == _SURFACE):
self.printOut("Data-Set has wrong type (UFO-Surface)", self.log.ERROR_TYPE, False)
return None
Imax = 0.
Cmax = 0
Rmax = 0.
headerV = self.Header.v
for c in range(0, self.points) :
d = self.Data[c]
r = np.asarray([d.v.x, d.v.y, d.v.z])
R = np.sqrt(np.dot(r-c, r-c))
if R > Rmax :
Rmax = R
if d.I > Imax:
Imax = d.I
Cmax = c
d = self.Data[Cmax]
c = np.asarray([d.v.x, d.v.y, d.v.z])
E = np.asarray([0.*1j, 0.*1j, 0.*1j])
for c in range(0, self.points) :
d = self.Data[c]
r = np.asarray([d.v.x, d.v.y, d.v.z])
R = np.sqrt(np.dot(r-c, r-c))
if R < Rmax*20. :
e = np.asarray([d.Er.x + d.Ei.x*1j, d.Er.y + d.Ei.y*1j, d.Er.z + d.Ei.z*1j])
E += np.sqrt(np.multiply(e, e.conj()).real )
pv = np.sqrt(np.multiply(E, E.conj()) ).real
pvn = np.sqrt(np.dot(pv, pv))
return pv/pvn
[docs]
def surfaceError ( self, rms = 1. ) :
"""
Implements a surface error with the rms given to the actual ufo surface.
Assumes a undisturbed perfectsurface data set (running the procedure twice results in sqrt(2) times rms.
Parameters:
rms (float): rms-value [mm] of the surface after the procedure
"""
self.printOut("Data-set %s : adding a surface RMS of %.3f mm" % (self.DataID, rms), self.log.FLOW_TYPE)
# loop all data points
for idx in range(0, self.points) :
displace = np.random.normal(scale = rms)
self.Data[idx].v.x += self.Data[idx].n.x * displace # and apply it to position vector
self.Data[idx].v.y += self.Data[idx].n.y * displace
self.Data[idx].v.z += self.Data[idx].n.z * displace
[docs]
def simulate ( self, target, wavelength, er = 1., tanD = 0., erOut = -1., trans = False):
"""
EM-simulation from the UFO data class (must be _CURRENT) to the target object (must be type _SURFACE).
The code will call the c-library with the code of the calculation to calculate the surface currents on the target object, resulting from the surface current induced E-Field at the location of the target.
Parameters:
target (UFODataClass): of type _SURFACE on which the current is calculated
wavelength (float): wavelength used for the simulation
er (float): permitivity of the medium in between the surfaces (is set automatically if data are available) -- no override
tanD (float): tangens delta (loss-factor) of the material in between surfaces (is set automatically if data are available) -- no override
erOut (float): permitivity of the medium after the target object; if negative, object is taken as a metal object (is set automatically if data are available) -- no override
Returns:
new (UFODataClass): with the simulation results
"""
self.printOut("Simulating e-field propagation via UFO c-code library:", self.log.FLOW_TYPE)
if (self.Wavelength != None) and wavelength < 0.:
wavelength = self.Wavelength
# create arrays
if self.Type != _CURRENT :
self.printOut("this class instance has no current information available => simulation not possible", self.log.ERROR_TYPE)
return None
new = deepcopy(target)
# set excitation current if not already available
if self.I0 == None :
self.setI0()
if not isinstance(self.er, (float, int, complex) ) :
self.er = er
if not isinstance(target.er, (float, int, complex) ) :
if erOut < 0 :
target.er = self.er
else :
target.er = erOut
if not isinstance(self.tanD, (float, int, complex) ) :
self.tanD = tanD
self.printOut("Wavelength : %.2f" % wavelength, self.log.DATA_TYPE, False)
self.printOut("Permittivity : %.2f" % self.er, self.log.DATA_TYPE, False)
self.printOut("Tanges-Delta : %.6f" % self.tanD, self.log.DATA_TYPE, False)
self.printOut("Permitivity after target : %.2f" % target.er, self.log.DATA_TYPE, False)
self.printOut("from", self.log.DATA_TYPE, False)
self.printOut("object : %s" % self.Name, self.log.DATA_TYPE, False)
v = [self.Header.v.x, self.Header.v.y, self.Header.v.z]
self.printOut("center : (%.2fmm, %.2fmm, %.2fmm)" %(v[0], v[1], v[2]), self.log.DATA_TYPE, False)
self.printOut("points : %d" % self.points, self.log.DATA_TYPE, False)
v = [target.Header.v.x, target.Header.v.y, target.Header.v.z]
self.printOut("to", self.log.DATA_TYPE, False)
self.printOut("object : %s" % target.Name, self.log.DATA_TYPE, False)
self.printOut("center : (%.2fmm, %.2fmm, %.2fmm)" %(v[0], v[1], v[2]), self.log.DATA_TYPE, False)
self.printOut("points : %d" % new.points, self.log.DATA_TYPE, False)
# simulate mirror
_cCalc_setLog( self.log )
#if self.SurfaceType == "Lens" : trans = True
new = _cCalc_CurrentDistFull(self, new, wavelength = wavelength, add = True, trans = trans )
new.Type = _CURRENT
new.Wavelength = wavelength
new.I0 = self.I0 # transfer initial excitation current troughout the simulations
new.er = target.er
# make sure field data are correct after lens simulation, assuming interface to vacuum, and skip impedance variations on angle of incident -- needs to be corrected!
new._modField ( np.sqrt( er ), onlyE = True )
if new.SurfaceType == "Lens" :
pass
return new