#!/bin/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/.
"""
# Package to simulate PAF-receiver systems
#
# The Package makes use of the UFO-package for simulating mirrors
#
#############################################################################
# import libs
#############################################################################
import sys
import glob
import subprocess
from copy import copy, deepcopy
import numpy as np
import matplotlib.pyplot as plt
from pathlib import Path
#############################################################################
# common functions
#############################################################################
_CLASSNAME = "UFO-Sim"
import LogClass
from Simulation import getColors
from Simulation import FarFieldDataClass
from Simulation.constants import _Tiso
#############################################################################
# common variables
#############################################################################
# initialize output class
log = LogClass.LogClass(_CLASSNAME)
#############################################################################
# shortcut to main classes (all making use of doc and simData)
#############################################################################
from .UFODataClass import UFOData
from . import UFO_CCode
from . import UFOGaussCalc
#############################################################################
#General methods
#############################################################################
[docs]
def setLog ( newLog ) :
global log
log = deepcopy( newLog )
log.addLabel ( cls = "UFO-Main" )
#####################################
# rotation of coordinate systems
[docs]
def getAngles ( zOrientation = [ 0., 0., 1.], xOrientation = [ 1., 0., 0. ] ) :
"""
Calculate rotation angles to get from a coordinate system with zAxis and xAxis given to the underlying coordinate system with
xAxis = [0, 0, 1], yAis = [0, 1, 0], and zAxis = [0, 0, 1]
Parameters:
zOrientation (array of float) : z-Axis directional vector
xOrientation (array of float) : x-Axis directional vector
Returns :
tuple (3 x 3 array of float) : rotation matrix into standart coordinate system
tuple (3 x 3 array of float) : rotation matrix from standart coordinate system to given coordinate system (inverse of the first return array)
tuple (3 element array of float) : rotation angles into standart coordinate system in degree [ rotation around x-axis, around y-axis, around z-axis ];
first rotating around x-, than around y- and last around z-axis corresponds to the first returned rotation matrix
"""
# calculate rotation angles to rotate z-axis vector to [0, 0, 1]
aX = np.arctan2( zOrientation[1], zOrientation[2] )
aY = -np.arctan2( zOrientation[0], np.sqrt( zOrientation[1]**2 + zOrientation[2]**2 ) )
# calculate x-vector orientation on x-y plane after de-rotating z-axis
irx = np.asarray([
[1, 0, 0],
[0, np.cos(aX), -np.sin(aX)],
[0, np.sin(aX), np.cos(aX)]
])
iry = np.asarray([
[ np.cos(aY), 0, np.sin(aY)],
[ 0, 1, 0],
[-np.sin(aY), 0, np.cos(aY)]
])
newX = np.dot( iry, np.dot( irx, xOrientation ) )
aZ = -np.arctan2(newX[1], newX[0])
irz = np.asarray([
[ np.cos(aZ), -np.sin(aZ), 0],
[ np.sin(aZ), np.cos(aZ), 0],
[ 0, 0, 1]
])
# calculate rotation matrix to new orientation
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]
])
rot = np.dot( rx, np.dot( ry, rz ) )
iRot = np.dot( irz, np.dot( iry, irx ) )
return iRot, rot, np.asarray( [ 180. * aX / np.pi, 180. * aY / np.pi, 180. * aZ / np.pi ] )
[docs]
def getRotMatrixForAxis ( axis = [ 1., 0., 0. ], angle = 0. ) :
"""
Calculate rotation matrix for a rotation around an axis given by its directional vector and the rotation angle.
Parameters:
axis (array of float) : rotation-axis directional vector
angle (float) : rotation angle [deg]
Returns:
tuple (3 x 3 array of float) : rotation matrix
"""
# normalize axis
axis = np.asarray( axis )
n = axis / ( np.sqrt( np.dot( axis, axis ) ) )
s = np.sin( np.pi * angle / 180. )
c = np.cos( np.pi * angle / 180. )
rot = np.asarray( [
[ n[0]**2 * ( 1 - c ) + c, n[0] * n[1] * ( 1 - c ) - n[2] * s, n[0] * n[2] * ( 1 - c ) + n[1] * s ],
[ n[1] * n[0] * ( 1 - c ) + n[2] * s, n[1]**2 * ( 1 - c ) + c, n[1] * n[2] * ( 1 - c ) - n[0] * s ],
[ n[2] * n[0] * ( 1 - c ) + n[1] * s, n[2] * n[1] * ( 1 - c ) + n[0] * s, n[2]**2 * ( 1 - c ) + c ]
] )
return rot
"""
Plot the geometrical layout including the mirror surfaces as three 2D projections.
Parameters:
objectList (list): of optical objects to be plotted
traceLine (numpy.array): with the ray path (e.g. output of GeometricTrace)
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)
cut (list): of 3-element lists. The first element in the sub-list denotes the axis (0 to 2) and the second element the boundary data, and the third the limit type (> 0 upper, ( < 0 lower, and == 0 no limit)
"""
log.printOut("Plot geometric 3D trace in projection including mirrors :", log.FLOW_TYPE)
trace = traceLine
if type(cut) == type([]) :
for c in cut :
trace = _cutTraceAxis( trace, c[1], c[0], c[2] )
colorList = getColors (len(objectList))
dataList = []
for obj in objectList :
dataList.append(obj.getData(True))
fig, (ax0, ax1, ax2) = plt.subplots(nrows=1, ncols=3, figsize=(12, 6))
ax0.set_ylabel("z-position [mm]")
ax0.set_xlabel("y-position [mm]")
ax0.set_aspect('equal', 'box')
for idx, data in enumerate(dataList) :
ax0.plot(data[:, 1], data[:, 2], color = colorList[idx], markersize = 1)
ax0.plot(trace [:, 1], trace [:, 2], "r-", label = "Beam-path" )
ax1.set_ylabel("z-position [mm]")
ax1.set_xlabel("x-position [mm]")
ax1.set_aspect('equal', 'box')
for idx, data in enumerate(dataList) :
ax1.plot(data[:, 0], data[:, 2], color = colorList[idx], markersize = 1)
ax1.plot(trace [:, 0], trace [:, 2], "r-")
ax2.set_ylabel("y-position [mm]")
ax2.set_xlabel("x-position [mm]")
ax2.set_aspect('equal', 'box')
for idx, data in enumerate(dataList) :
ax2.plot(data[:, 0], data[:, 1], color = colorList[idx], markersize = 1)
ax2.plot(trace [:, 0], trace [:, 1], "r-")
if figname != "" :
fig.savefig(figname)
if show :
plt.show()
else :
plt.show()
plt.close()
#####################################
# physical optics tracing
[docs]
def PhysicalTrace ( objects, wavelength, lossFactor = 0. , lossTemperature = _Tiso, waist = None, phCenter = [] ) :
"""
Physical simulation of a sequence of mirrors and return the last object with the current distribution and e-field.
Parameters:
objects (list): of UFOData-class instances in the order as they shall be traced. The first element must have a pre-defined current distribution
wavelength (float): wavelength of the optics setup (used to determine the opening angle of the initial beam)
lossFactor (float): 0 to 1., denoting the losses in the field during the trace.
lossTemperature (float): temperature [K] at which losses are dumped
waist (float or None): If given, a waist of this size is calculated on the first object of the objects-list
Returns:
tuple (tuple): curr, loss, lossT (UFO-Object with the result data (current type))
"""
log.printOut( "Physical optics trace at wavelength %.2f mm :" % wavelength, log.FLOW_TYPE)
if len(objects) < 2 :
log.printOut( " not enough objects given for a simulation, need in minimum two!", log.WARNING_TYPE, False )
return None, None, None
if isinstance ( objects[0], FarFieldDataClass.FarFieldData ) : # check if first element is far-field data
phCenter = np.asarray( [objects[1].Header.v.x, objects[1].Header.v.y, objects[1].Header.v.z ] ) - np.asarray( objects[0].Center )
objects[1] = objects[0].calculateCurrent( objects[1], add = False ) # simulate far-field pattern to first mirror
objects.pop(0) # remove first element (far-field pattern) from simulation list
else : # or UFO object
if isinstance(waist, ( float, int ) ) :
objects[0].makeWaist(waist)
phCenter = []
# calculate beam-direction
phCenter = np.asarray( [objects[1].Header.v.x, objects[1].Header.v.y, objects[1].Header.v.z ] ) - np.asarray( [objects[0].Header.v.x, objects[0].Header.v.y, objects[0].Header.v.z ] )
objects[0].Wavelength = wavelength
curr = objects[0]
loss, lossT = curr.applyLosses( lossFactor, lossTemperature, phCenter )
for idx in range(1, len(objects)):
curr = curr.simulate(objects[idx], wavelength)
return curr, loss, lossT
#####################################
# plotting
[docs]
def plotOpticalCircuit ( objects, startWaist, wavelength, startDistance = [], size = 4., figname = "", show = False, maxX = None ):
"""
Plot of the optical circuit of the objects in line.
Parameters:
objects (list): of UFOData-class instances in the order as they shall be traced. the first element is the transmitting surface
startWaist (float): or list of floats, waist size(s) of the starting waist(s)
wavelength (float): or list of floats, wavelength of the optics setup tracing
startDistance (float): or list of floats giving the difference of the start-waist distance from object[0] (optional)
size (float): denoting the multiplicator of the contour. e.g. a size of 5 results in the 5 x wz contour in diameter!
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)
maxX (boolean): maximum limit in x-direction (to cut of a plot)
"""
#################
# internal methods
def plotWaist(ax, waist, wl, z_start, z_stop, x0, size, color = "black") :
"""
Internal method to plot the waist contour on axis ax.
Parameters:
ax (axis object ?): for the plotting
waist (float): waist size [mm]
wl (float): wavelength [mm]
z_start (float): start-point in waist coordinate system (waist is at z == 0)
z_stop (float): stop-point in waist coordinate system (waist is at z == 0)
x0 (float?): absolute coordinate offset for plotting
size (float?): contour to plot (1: 1/e contour, 2: 2 time 1/e contour, ...)
color (string): plot-color
"""
z_step = 1.
if z_start > z_stop : _z_step = -1.
contour = []
x = []
actX = x0
for z in np.arange(z_start, z_stop, z_step) :
zr = np.pi*waist**2./wl
contour.append(size*waist*np.sqrt(1+z**2./zr**2)/2.)
x.append(actX)
actX += 1.
contour = np.asarray(contour)
line, = ax.plot(x, contour, color = color)
line, = ax.plot(x, -contour, color = color)
return actX, line
# end of method plotWaist
# end of internal method definitions
#################
log.printOut("Plot the optical circutry :", log.FLOW_TYPE)
# transform input data formats
if isinstance(startWaist, (float, int)) :
startWaist = [startWaist]
if isinstance(wavelength, (float, int) ) :
wavelength = [wavelength]
if isinstance(startDistance, (float, int) ) :
startDistance = [startDistance]
if len(startDistance) == 0 :
startDistance = [0.]*len(startWaist)
if ( len(wavelength) != len(startWaist) ) or ( len(startWaist) != len(startDistance) ) :
return
log.printOut(" using %d objects" % len(objects), log.DATA_TYPE, False)
log.printOut(" and %d wavelength" % len(wavelength), log.DATA_TYPE, False)
# establish color-list
colors = ['green', 'blue', 'magenta', 'purple', 'red', 'orange', 'yellow']
while len(colors) < len(wavelength) : colors += colors
# extract mirror distances
start = False
distances = []
for obj in objects :
v = np.asarray([obj.Header.v.x, obj.Header.v.y, obj.Header.v.z])
if not isinstance(start, bool) :
distances.append(np.sqrt(np.dot(v-start, v-start)))
start = v
ledgendLabels = [] # ledgend label
ledgendLines = [] # and ledgend line objects
# open plot-screen
fig, ax = plt.subplots(figsize=(16, 10))
ax.set_ylabel ("Beam-size (%.1f $\\times\\omega_0$) and mirror size [mm]" % (size/2.))
ax.set_xlabel ("Distance [mm] from nominal horn Waist at $\\lambda = %.1f mm$" % wavelength[0])
distances_org = deepcopy(distances)
for w_idx, wl in enumerate(wavelength) : # loop all wavelength given
log.printOut(" - Using wavelength %.2fmm (index = %d)" % (wl, w_idx), log.DEBUG_TYPE, False)
distances = deepcopy(distances_org)
distances[0] = distances[0] - startDistance[w_idx]
mw1 = startWaist[w_idx] # start waist
md1 = distances[0] # start input waist distance
cstart = 0. # contour plotting start value (z-coordinate in the waist system)
cstop = md1 # contour plotting stop value (z-coordinate in the waist system)
xAct = startDistance[w_idx] # actual x-position in the plot
posList = []
maxHeight = 0.
# loop all objects and create the ray-tracing
for idx, obj in enumerate(objects) :
# plot beam-contour
if idx > 0 : # first object represents the horn antenna (or at least a surface with the start-waist current contribution)
if obj.FocalLength != None : # if a focal length is given, calculate mirror characteristics for the actual wavelength
mw1, mw2, md1, md2, mf = UFOGaussCalc.complete( wl, w1 = mw1, w2 = None, d1 = md1, d2 = None, f = obj.FocalLength)
cstart = -md2 # adapt start point for the contour
mw1 = mw2 # make output waist being the input waist for the next mirror
if idx < len(objects) -1 : # if it is not the last mirror
md1 = distances[idx] - md2 # calculate input waist distance for the next active mirror
else : # if no focal length given, treate the surface as a flat mirror
if idx < len(objects) -1 : # if it is not the last object
md1 += distances[idx] # add following distance to input length of the next active mirror
if idx == len(objects) - 1: # if it is the last object
cstop = 0. # stop contour at the output waist of the last mirror
if isinstance(maxX, float) :
if maxX > xAct :
cstop = cstart + maxX-xAct
else:
cstop = cstart + distances[idx] # else stop contour at the next mirror
# plot objects
if w_idx == 0 :
log.printOut(" plotting object %d with surface type %s" % (idx, obj.SurfaceType), log.DEBUG_TYPE, False)
# get x-dimension for the plot-size
o = deepcopy(obj)
o.resetPosition()
xmax = 0.
for d in o.Data :
if d.v.x**2. > xmax : xmax = d.v.x**2.
xmax = np.sqrt(xmax)
radX = []
radY = []
posList.append(xAct)
if xmax > maxHeight :
maxHeight = xmax
color = "black"
if (obj.SurfaceType == "Mirror") or ( (obj.SurfaceType == "Lens") and ( obj.Name.find("(2)") < 0 ) ) : # is object an active mirror or lens (Only plot first surface of the lens)?
# if ( obj.Name.find("(2)") < 0 ) : # second lens surface should not be plotted
# if isinstance(obj.FocalLength, float) or isinstance(obj.FocalLength, int): # second lens surface has no focal length given and should not being plotted
xv = xmax
if xv > obj.FocalLength*2. : xv = obj.FocalLength*2.
angle = np.arcsin(xv /obj.FocalLength/2.)
for a in np.arange(-angle, angle, angle/30.) :
radX.append(np.cos(a)*2.*obj.FocalLength)
radY.append(np.sin(a)*2.*obj.FocalLength)
if len(radX) == 1 :
ref = radX[0]
radX[-1] -= ref
color = "red"
elif ( obj.SurfaceType == "Flat" ) or ( obj.SurfaceType == "Aperture" ): # is object a flat mirror or an aperture?
radX = [0., 0.]
radY = [-xmax, xmax]
color = "blue"
log.printOut(" surface points x : %.2f, %.2f" % (radY[0], radY[1] ), log.DATA_TYPE, False)
if (len(radX) > 0) and (idx > 0) :
radXo = np.asarray(radX)
radYo = np.asarray(radY)
radXo = np.concatenate((radXo, -1.*radXo))
radYo = np.concatenate((radYo, 1.*radYo))
RefAngle = obj.RefAngle
if ( RefAngle == None ) or ( RefAngle == 180. ) :
RefAngle = 0.
RefAngle = -np.pi * RefAngle / 360.
radX = radXo*np.cos( RefAngle ) - radYo*np.sin( RefAngle )
radY = radXo*np.sin( RefAngle ) + radYo*np.cos( RefAngle )
ax.plot( radX + xAct, radY, color = color)
ax.plot([xAct, xAct],[radY[0], -radY[0]+20], linestyle = "dashed", color = color)
NameString = obj.Name
bracketString = ""
if (obj.FocalLength != None ):
bracketString += " f = %.1f mm" % obj.FocalLength
addedF = True
if ( obj.RefAngle != None ) and ( obj.RefAngle != 0. ):
if bracketString != "" :
bracketString += ",\n "
bracketString += " reflection angle = %.1f deg." % obj.RefAngle
if bracketString != "" :
bracketString = " ("+ bracketString + ")"
ax.text(xAct, -radY[0]+20, NameString + bracketString, rotation=45, color = color)
if idx == 0 :
startWaist[w_idx]
hornx = (np.asarray([0, 0, -350, -350, -370, -370]) + 50.) /13.5*wl
horny = np.asarray([0., size/2.*startWaist[w_idx], 10/13.5*wl, 7/13.5*wl, 7/13.5*wl, 0])
if np.amin(hornx) > -150 :
posText = np.amax(horny)
else :
posText = 0
ax.plot(hornx, horny*1.2, "k-")
ax.plot(hornx, -horny*1.2, "k-")
ax.text( -80, posText, "Horn\nantenna", color = "black", ha = "right")
# plot contour and update actual x-position
xAct, line = plotWaist(ax, mw1, wl, cstart, cstop, xAct, size, color = colors[w_idx])
# reset start point of the contour to the actual countour end
cstart = cstop
# end of object loop
for idx, pos in enumerate(posList) :
ax.plot([pos,pos], [0., -maxHeight-50], linestyle = "dashed", color = "black")
if idx < len(posList)-1 :
if distances[idx] > 0. : # take care of second lens-surface
ax. arrow(pos + distances[idx]/2., -maxHeight-30, -distances[idx]/2.+10., 0., color = "black", width = 2)
ax. arrow(pos + distances[idx]/2., -maxHeight-30, distances[idx]/2.-10., 0., color = "black", width = 2)
ax.text(pos + distances[idx]/2., -maxHeight-50, "%.1f mm" % distances[idx], ha = 'center')
# add ledgend
ledgendLabels.append("$%.1f \\times \\omega_z(z) for \\lambda = %.1f mm$" %(size/2., wl))
ledgendLines.append(line)
#end of wavelength loop
xmin, xmax = ax.get_xlim()
xmax *= 1.2
if isinstance(maxX, float) : xmax = maxX
ax.set_xlim([xmin, xmax])
ymin, ymax = ax.get_ylim()
ax.set_ylim([ymin, ymax*1.75])
ax.legend(ledgendLines, ledgendLabels, loc = 'upper left' )
if figname != "" :
fig.savefig(figname)
if show :
plt.show()
else :
plt.show()
plt.close()