#!/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/.
"""
#############################################################################
# import libs
#############################################################################
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.mlab as ml
from copy import copy, deepcopy
from pathlib import Path
try:
from scipy.interpolate import griddata
_SciPy = True
except:
_SciPy = False
import LogClass
from Simulation.UFO import UFODataClass
from Simulation import CSTDataClass
from Simulation import makeMovie
from Simulation.constants import _ETA, _C, _I0iso, _Bandwidth, _Impedance
from Simulation.constants import _GridSize
from Simulation import getColors
#############################################################################
# constants
#############################################################################
#
_CLASSNAME = "FarFieldData"
#############################################################################
# classes
#############################################################################
[docs]
class FarFieldData :
"""
Class to handle far field data
"""
#####################################
# internal data and init method
# Keywords for the file
__KEYS = {
"Type" : 'FARFIELD:' ,
"Wavelength" : 'WAVELENGTH:' ,
"Calculated" : 'CALCULATED:' ,
"Comment" : '#' ,
"Center" : 'CENTER:' ,
"Ox" : 'ORIENTATION-X:' ,
"Oz" : 'ORIENTATION-Z:' ,
"PhaseCenter" : 'PHASECENTER:' ,
"Gain" : 'GAIN:' ,
"Phase" : 'PHASE:' ,
"Data" : 'DATA:' ,
"ReferenceDist" : 'REFERENCE-DIST:' ,
"I0" : 'I0:' ,
"Bandwidth" : 'Bandwidth:' ,
"S11" : 'S11:' ,
"S21" : 'S21:' ,
"Impedance" : 'IMPEDANCE:'
}
__Type = {
"NA" : 'Unknown' ,
"GAUSS" : 'Gaussian' ,
"DIPOLE" : 'Dipole' ,
"CST" : 'CST'
}
# init method
def __init__ ( self, wavelength, raster = 5, log = None ):
"""
Initialize class. A wavelength can be given.
"""
self.setLog ( log )
self.TypeDict = self.__Type
# Header data
self.Wavelength = wavelength # Wavelength in mm
self.Type = "" # Far-Field data type
# Koordinate system definitions (used to rotate and shift the data in 3D space
self.OrientationX = [1., 0., 1.]
self.OrientationZ = [0., 0., 1.] # direction of propagation and direction with Theta = 0
self.Center = [0., 0., 0.] # Center of the Far-field data in 3D space [mm]
self.PhaseCenter = [0., 0., 0.] # Center of the phase-distribution, relative to the coordinate center [mm]
self.ReferenceDist = 1000. # default reference distance of the far-field pattern [mm]
self.Rot = np.identity(3) # Rotation matrix for the output data
self.invRot = np.linalg.inv(self.Rot) # inverse rotation
self.raster = raster # raster of the farfield in degrees
# Field data
self.Gain = 1. # Gain to boost the field if required
self.Phase = 0. # Phase of the file data [radian]
# Exitation data
self.I0 = 1. # excitation current [A]
self.Impedance = _Impedance # antenna impedance [Ohm]
self.Bandwidth = _Bandwidth # calculation bandwidth of the far-field pattern
# S-Parameters
self.S11 = 0. # S11
self.S21 = 1. # S21
# data array
self.FarField = [] # 2D List of far-field data with index [index-Phi][index-Theta]
# each Array entry contains a List with Phi, Theta, E_Phi, and E_Theta values
self.freqDict = {} # dictionary with frequencies and belonging file-names
self.cst = CSTDataClass.CST( log = self.log ) # cst class instance allowing to access CST simulation data directly
[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, noSpace = True, lineFeed = True ) :
"""
Wrapper for the printOut method of the actual log class.
"""
self.log.printOut(line, message_type, noSpace, lineFeed)
#####################################
# operator overloading
def __add__ ( self, FFClass ) : # overload adding far-field data with combining them
if (len(self.FarField) < 2) or (len(FFClass.FarField) < 2) : # check if data are available, if not return empty far-field
self.printOut("Not enough data available, returning empty class", self.log.WARNING_TYPE)
return FarFieldData()
if (len(self.FarField[0]) < 2) or (len(FFClass.FarField[0]) < 2) : # check if data are available, if not return empty far-field
self.printOut("Not enough data available, returning empty class", self.log.WARNING_TYPE)
return FarFieldData()
dphi = np.absolute(self.FarField[1][0][0] - self.FarField[0][0][0])
dtheta = np.absolute(self.FarField[0][1][1] - self.FarField[0][0][1])
return self.combineFarField(FFClass, dphi, dtheta)
def __sub__ ( self, FFClass ) : # overload substracting far-field data by shifting phase of the FFClass by 180
FFClassMinus = deepcopy(FFClass)
FFClassMinus.setPhase(180.)
return self.__add__(FFClassMinus)
#####################################
# File-IO methods
[docs]
def writeFile ( self, filename ):
"""
Method to write the Far-Field data to a data file.
Parameters:
filename (string): Name of the file to be written.
"""
self.printOut("Writing File %s" % filename, self.log.FLOW_TYPE)
fh = open(filename, 'w') # open file for writing
# write header
fh.write(self.__KEYS['Comment'] + " UFO FarField file\n" )
fh.write(self.__KEYS['Comment'] + " \n" )
fh.write(self.__KEYS['Comment'] + " Data-Header\n" )
fh.write(self.__KEYS['Type'] + " # Type Identifier (only lines after Type-ID are parsed during reading\n" )
fh.write(self.__KEYS['Wavelength'] + " %.3f mm\n" % self.Wavelength )
fh.write(self.__KEYS['Calculated'] + " %s\n" % self.Type )
fh.write(self.__KEYS['Comment'] + " \n" )
fh.write(self.__KEYS['Comment'] + " Koordinate definitions:\n" )
fh.write(self.__KEYS['Center'] + " %.5f %.5f %.5f\n" % (self.Center [0], self.Center [1], self.Center [2]) )
fh.write(self.__KEYS['Ox'] + " %.5f %.5f %.5f\n" % (self.OrientationX[0], self.OrientationX[1], self.OrientationX[2]) )
fh.write(self.__KEYS['Oz'] + " %.5f %.5f %.5f\n" % (self.OrientationZ[0], self.OrientationZ[1], self.OrientationZ[2]) )
fh.write(self.__KEYS['PhaseCenter'] + " %.5f %.5f %.5f\n" % (self.PhaseCenter [0], self.PhaseCenter [1], self.PhaseCenter [2]) )
fh.write(self.__KEYS['ReferenceDist'] + " %.5f\n" % self.ReferenceDist )
fh.write(self.__KEYS['Comment'] + " \n" )
fh.write(self.__KEYS['Comment'] + " Signal modifications:\n" )
fh.write(self.__KEYS['Gain'] + " %.5f\n" % self.Gain )
fh.write(self.__KEYS['Phase'] + " %.5f deg\n" % self.Phase )
fh.write(self.__KEYS['Comment'] + " Exitation data (I0 and input impedance)\n")
fh.write(self.__KEYS['I0'] + " %.5f A\n" % self.I0)
fh.write(self.__KEYS['Bandwidth'] + " %.2f Hz\n" % self.Bandwidth)
fh.write(self.__KEYS['S11'] + " %.5f\n" % self.S11)
fh.write(self.__KEYS['S21'] + " %.5f\n" % self.S21)
fh.write(self.__KEYS['Impedance'] + " %.5f Ohm\n" % self.Impedance)
fh.write(self.__KEYS['Comment'] + " \n" )
fh.write(self.__KEYS['Comment'] + " Data-Block (index-Phi index-Theta Phi[deg] Theta[deg] E-Phi_real E-Phi_imag E-Theta_real E-Theta_imag):\n" )
for indexPhi in range(0, len(self.FarField)) :
F = self.FarField[indexPhi]
for indexTheta in range(0, len(F)) :
fh.write(self.__KEYS['Data'] + " %d %d %.6f %.6f %.10f %.10f %.10f %.10f\n"
% ( indexPhi, indexTheta, F[indexTheta][0], F[indexTheta][1],
F[indexTheta][2].real, F[indexTheta][2].imag, F[indexTheta][3].real, F[indexTheta][3].imag) )
fh.write(self.__KEYS['Comment'] + " \n" )
fh.write(self.__KEYS['Comment'] + " end of file\n" )
fh.close()
[docs]
def readUFOFarfieldFile ( self, filename ):
"""
A method to read a farfield file in UFO (?) format.
This fileformat seems to correspond to the one found in the writeFile-Method of this class.
Parameters:
filename (string): the path of the file to be read
Returns:
idOK (boolean): whether or not reading the file was successful
"""
self.printOut("Reading Farfield-file in UFO format %s" % filename, self.log.FLOW_TYPE)
idOK = False
with open(filename, 'r') as f : # open file for reading
for line in f : # and loop all lines in the file
line = line.strip() # remove CR / NL / EOF from line end
columns = line.split() # split line
if len(columns) > 0:
if columns[0] == self.__KEYS.get('Type'): # check for correct type-ID
idOK = True
elif columns[0] == self.__KEYS.get('Wavelength') :
self.Wavelength = float(columns[1])
elif columns[0] == self.__KEYS.get('Calculated') :
self.Type = columns[1]
elif columns[0] == self.__KEYS.get('Center') :
self.Center[0] = float(columns[1])
self.Center[1] = float(columns[2])
self.Center[2] = float(columns[3])
elif columns[0] == self.__KEYS.get('Ox') :
self.OrientationX[0] = float(columns[1])
self.OrientationX[1] = float(columns[2])
self.OrientationX[2] = float(columns[3])
elif columns[0] == self.__KEYS.get('Oz') :
self.OrientationZ[0] = float(columns[1])
self.OrientationZ[1] = float(columns[2])
self.OrientationZ[2] = float(columns[3])
elif columns[0] == self.__KEYS.get('PhaseCenter') :
self.PhaseCenter[0] = float(columns[1])
self.PhaseCenter[1] = float(columns[2])
self.PhaseCenter[2] = float(columns[3])
elif columns[0] == self.__KEYS.get('ReferenceDist') :
self.ReferenceDist = float(columns[1])
elif columns[0] == self.__KEYS.get('Gain') :
self.Gain = float(columns[1])
elif columns[0] == self.__KEYS.get('Phase') :
self.Phase = float(columns[1])
elif columns[0] == self.__KEYS.get('I0') :
self.I0 = float(columns[1])
elif columns[0] == self.__KEYS.get('Bandwidth') :
self.Bandwidth = float(columns[1])
elif columns[0] == self.__KEYS.get('S11') :
self.S11 = float(columns[1])
elif columns[0] == self.__KEYS.get('S21') :
self.S21 = float(columns[1])
elif columns[0] == self.__KEYS.get('Impedance') :
self.Impedance = float(columns[1])
elif columns[0] == self.__KEYS.get('Data') :
indexPhi = int(float(columns[1]))
indexTheta = int(float(columns[2]))
Phi = float(columns[3])
Theta = float(columns[4])
E_Phi = float(columns[5]) + 1j*float(columns[6])
E_Theta = float(columns[7]) + 1j*float(columns[8])
while len(self.FarField) < indexPhi+1 :
self.FarField.append([])
while len(self.FarField[indexPhi]) < indexTheta+1 :
self.FarField[indexPhi].append([])
self.FarField[indexPhi][indexTheta] = [Phi, Theta, E_Phi, E_Theta]
f.close()
return idOK
[docs]
def readCSTFarfieldFile ( self, filename, cstImpedance = _Impedance ) :
"""
Method to read a CST ascii file with Far-Field data.
Parameters:
filename (string): Name of the input file
single (boolean?): Handling a folder setup (False) or read a single file (True)
Returns:
idOK (boolean): error flag (true on success, false on error)
"""
self.printOut("Reading Farfield-file in CST format %s" % filename, self.log.FLOW_TYPE)
dataValid = False # data valid flag
lcount = 0 # line-counter
lExp = -10 # line with experimental data (Power, frequency)
lNum = -10 # line with number of elements data
lCenter = -10 # line with center position
lzAxis = -10 # line with z-axis data
lxAxis = -10 # line with x-axis data
nPhi = 0 # number of Phi values
nTh = 0 # number of Theta values
with open(filename, 'r') as f : # open file
for line in f : # and loop all lines in the file
line = line.strip() # remove CR / NL / EOF from line end
if not dataValid : # if line is not a data-line parse it for key-words
if line == "// Radiated/Accepted/Stimulated Power , Frequency" : # parse line contend for experimental data
lExp = lcount # set line-number of the frequency entry
if line == "// Position" : # parse line contend for center position
lCenter = lcount+1 # set line-number of the center position
if line == "// xAxis" : # parse line contend for center position
lxAxis = lcount+1 # set line-number of the center position
if line == "// zAxis" : # parse line contend for center position
lzAxis = lcount+1 # set line-number of the center position
if line == "// >> Total #phi samples, total #theta samples" : # parse line contend for number of Phi & Theta values
lNum = lcount + 1 # set line number of the sampling data
if line == "// >> Phi, Theta, Re(E_Theta), Im(E_Theta), Re(E_Phi), Im(E_Phi):" : # parse line contend for start of the data block
dataValid = True # indicate start of the data-block
for i in range(nPhi): # pre-define 2d-array
row = []
for j in range(nTh):
row.append([])
self.FarField.append(row)
else : # otherwise data are valid
columns = line.split() # extract values
c0 = float(columns[0])
c1 = float(columns[1])
iPh = int ( (nPhi-1)*float(c0)/360.) # calculate array index for Phi and Theta
iTh = int ( (nTh -1)*float(c1)/180.)
self.FarField[iPh][iTh] = [ c0, c1, # compose complex Ephi and Eth-data
(float(columns[4])+1j*float(columns[5])) , # (Eth comes first in the data-line, but array
(float(columns[2])+1j*float(columns[3])) ] # is [Phi, Theta, Ephi, Eth]!)
# extract experimental data
if (lcount > lExp) and (lcount < lExp+5) : # at correct line-number extract experimental data
columns = line.split()
if lcount - lExp == 1 :
self.S21 = float(columns[0]) # extract S21
if lcount - lExp == 2 :
self.S11 = float(columns[0]) # extract S11
if lcount - lExp == 3 :
self.I0 = np.sqrt(2.*float(columns[0])/cstImpedance) # extract I0
self.S21 = np.sqrt(self.S21 / float(columns[0]))
self.S11 = np.sqrt(1. - self.S11 / float(columns[0]))
if lcount - lExp == 4 :
self.Wavelength = _C / float(columns[0]) * 1000. # or Wavelength in mm
# extract number of data-points
if lcount == lNum : # at correct line-number extract array dimensions
columns = line.split()
nPhi = int(columns[0])
nTh = int(columns[1])
# extract center position
if lcount == lCenter : # at correct line-number extract center positon
columns = line.split()
self.PhaseCenter = [ float(columns[0])*1000., float(columns[1])*1000., float(columns[2])*1000. ]
# extract x-axis orientation
if lcount == lxAxis : # at correct line-number extract x-axis
columns = line.split()
self.OrientationX = [ float(columns[0]), float(columns[1]), float(columns[2]) ]
# extract z-axis orientation
if lcount == lzAxis : # at correct line-number extract z-axis
columns = line.split()
self.OrientationZ = [ float(columns[0]), float(columns[1]), float(columns[2]) ]
lcount += 1 # increment line-counter
f.close()
return dataValid
[docs]
def importFarfieldFile ( self, filename, setPhaseCenter = False, single = True, fileformat = "" ):
"""
Method to import a Farfield data file in either CST or UFO format
Parameters:
filename (string): Name of the input file
setPhaseCenter (boolean): Use phase center (or fit it in case of CST-files) -- True; or set it to zero
single (boolean): Handling a folder setup (False) or read a single file (True)
fileformat (string): File-Format indicator, empty string "": try out, '"UFO": internal fime format, "CST": CST fileformat
Returns:
idOK (boolean): error flag (true on success, false on error)
"""
self.FarField = [] # delete all data
idOK = False
if fileformat in ( "UFO", "" ):
idOK = self.readUFOFarfieldFile(filename)
if not setPhaseCenter :
self.PhaseCenter = [0, 0, 0]
if ( fileformat == "CST" ) or ( ( fileformat == "" ) and ( idOK == False ) ):
idOK = self.readCSTFarfieldFile(filename)
self.moveCenter (0., 0., 0.) # absolute center not given, move it to zero
self.movePhaseCenter(0., 0., 0.)
if setPhaseCenter :
self.estimatePhaseCenter() # estimate phase-center from far-field distribution
self.ReferenceDist = 1000. # reference distance in CST is 1m
self.Gain = 1.
self.Phase = 0.
self.Bandwidth = 1000000. # default CST calculation bandwidth (1MHz)
self.Type = self.__Type['CST'] # set calculation type
self.adjustFarfield()
# For a folder setup the frequency dictionary was set with calling the "setFolder" method.
# When reading a single file, the frequency dictionary must be set correctly
if single :
# delete key-list
self.freqDict = {}
frequency = _C/self.Wavelength/1.e6 # frequency [GHz]
self.freqDict["%.2f" % frequency] = filename
if not idOK:
self.printOut("Error while reading file", self.log.ERROR_TYPE, False)
else : # adapt excitation current to the nominal current of the simulator
self.scaleI0()
return idOK
[docs]
def adjustFarfield ( self ):
"""
Re-calculate internal rotation matrices, used after importing data including coordinate system (implemented into import function above).
"""
# get rotation angles to get back to positive Z-axis being aligned with theta = 0
alpha_z = - 180./np.pi * np.arctan2( self.OrientationZ[1], self.OrientationZ[0])
alpha_y = - 180./np.pi * np.arctan2( np.sqrt( self.OrientationZ[0]**2 + self.OrientationZ[1]**2), self.OrientationZ[2] )
# and de-rotate the fields
self.rotate(0., 0., alpha_z)
self.rotate(0., alpha_y, 0. )
# check x-axis angle to get rotation around surface normal vector
alpha_n = 180./np.pi * np.arctan2(self.OrientationX[1], self.OrientationX[0])
# and de-rotate the fields
self.rotate(0., 0., alpha_n )
# reset rotation matrices and orientation vectors
self.Rot = np.identity(3) # Rotation matrix for the output data
self.invRot = np.linalg.inv(self.Rot) # inverse rotation
self.OrientationX = [1., 0., 0.]
self.OrientationZ = [0., 0., 1.]
# re-calculate matrices and vectors from rotation angles determined above
self.rotate(0., 0., -alpha_n)
self.rotate(0., -alpha_y, 0. )
self.rotate(0., 0., -alpha_z)
[docs]
def setCST ( self, CSTproject, folder = "", fileNameBase = "ffx_", PolAngle = 0. ) :
"""
Read CST-files from the project, store them in an internal format to folder 'folder' and select this folder for data access.
Parameters:
CSTproject (string): name of the CST project to access
folder (string): name of the folder to save the far-field sources (if empty, no access to far-field source data)
fileNameBase (string): base-name for the new farfield files. Available farfields will be numbered in sequence.
PolAngle (float): gives th polarization angle [deg] of the created source (0 deg = along x-axis)
"""
self.printOut("Setting access to CST project '%s'" % CSTproject, self.log.FLOW_TYPE)
self.cst = CSTDataClass.CST( CSTproject, log = self.log )
if folder != "" :
self.printOut("trying to access the far-field-source files", self.log.DEBUG_TYPE, False)
ok, ffs = self.cst.readFarFieldSources()
if ok:
for i, ff in enumerate(ffs) :
ff.rotate(0, 0, PolAngle)
ff.writeFile(folder + "/" + "%s%d.ffs" % (fileNameBase, i))
self.setFolder(folder, wildcard = fileNameBase + "*.ffs")
else:
self.printOut("No far-field data available!", self.log.WARNING_TYPE, False)
[docs]
def setFolder ( self, folder, wildcard = "*.*" ) :
"""
Set a folder with far-field pattern for various frequencies.
Parameters:
folder (string?): name of the folder to be used as far-field source
wildcard (string?): wildcard of the filename
Returns:
list of frequencies available
"""
self.printOut("Setting folder with far-field data to '%s'" % folder, self.log.FLOW_TYPE)
self.printOut("using file-name wildcard '%s'" % wildcard, self.log.DATA_TYPE, False)
# delete key-list
self.freqDict = {}
# get all file-names in the folder
fileList = [f"{f}" for f in Path(folder).glob(wildcard)]
# loop all files and extract frequency
for fileName in fileList :
ok = self.importFarfieldFile(fileName, single = False)
if ok :
frequency = _C/self.Wavelength/1.e6 # frequency [GHz]
self.freqDict["%.2f" % frequency] = fileName
self.printOut("-- '%s' is a far-field file at frequency %.2fGHz (%.2f mm wavelength)" %(fileName, frequency, self.Wavelength), self.log.DEBUG_TYPE, False)
else :
self.printOut("-- '%s' is not a far-field file" % fileName, self.log.DEBUG_TYPE, False)
return sorted(self.freqDict)
[docs]
def setFrequency ( self, frequency ) :
"""
Loads the far-field file belonging to the given frequency (closest is chosen).
Parameters:
frequency (float): frequency [GHz] of the far-field pattern to be loaded (it is the nearest available frequency chosen!)
"""
self.printOut("Selecting frequency %.2fGHz" % frequency, self.log.DEBUG_TYPE)
if self.Type == self.__Type['DIPOLE'] :
self.printOut("changing frequency of the Dipole profile", self.log.DEBUG_TYPE, False)
self.Wavelength = _C / frequency /1.e6
self.createFarfield("Dipole", intensity = self.intensity, width = self.width, antennaLength = self.antennaLength)
elif self.Type == self.__Type['GAUSS'] :
self.printOut("changing frequency of the Gaussian profile", self.log.DEBUG_TYPE, False)
self.Wavelength = _C / frequency /1.e6
elif len(self.freqDict) > 0 :
self.printOut("using file from file list", self.log.DEBUG_TYPE, False)
dist = 99999999
select = 0
for f in sorted(self.freqDict) :
d = np.absolute(frequency - float(f))
if d < dist:
select = f
dist = d
self.importFarfieldFile(self.freqDict[select], single = False)
else :
self.printOut("no folder selected or no far-field data available in the folder", self.log.ERROR_TYPE)
#####################################
# Methods to establish a distribution
[docs]
def createWaist ( self, waist = 10, intensity = 1, stepsPerWaist = 50):
"""
creates a Gaussian far field distribution corresponding to the far field of the given waist size.
"""
width = 180 / np.pi * self.Wavelength / waist / np.pi
thetaMax = 5 * width
if thetaMax > 180 : thetaMax = 180
self.raster = width / stepsPerWaist
self.createFarfield(Type = "Gaussian", intensity = intensity, width = width, thetaMax = thetaMax)
[docs]
def calculateGaussianFarfieldPoint ( self, Phi, Theta, intensity = 1, width = 90 ) :
"""
Calculates one point (Phi,Theta) of a farfield for a Gaussian distribution.
With linear polarization vector in cross-Theta direction.
Parameters:
Phi (float): Phi angle of data point
Thete (float): Theta angle of data point
intensity (float): intensity of the distribution
width (float): 1/e width of the resulting field (in degrees)
Returns:
list (list): phi angle, theta angle, Phi component of farfield, Theta component of farfield
"""
gauss = intensity*np.exp(-(Theta/width)**2) # calculate intensity for Phi/Theta direction
if Theta > 90 : gauss = 0 # use only forward hemisphere
E_theta = np.sin(Phi/180*np.pi)*gauss # calculate e-field in polar coordinates
E_phi = np.cos(Phi/180*np.pi)*gauss
return [Phi, Theta, E_phi+0j, E_theta+0j] # append new data point
[docs]
def calculateDipoleFarfieldPoint ( self, Phi, Theta, antennaLength = 100) :
"""
Calculates one point (Phi,Theta) of a farfield for a dipole with groundplane.
The dipole is aligned with the X-axis and has a distance of wavelength/4 from the groundplane.
Parameters:
Phi (float): Phi angle of data point
Thete (float): Theta angle of data point
antennaLength (float) : length of the dipole section (oriented along the y-axis)
Returns:
list (list): phi angle, theta angle, Phi component of farfield, Theta component of farfield
"""
ratio = antennaLength / self.Wavelength
k0 = 2*np.pi/self.Wavelength
pol_vector = [1,0,0] # element orientation [x,y,z]
zOffset = antennaLength / 4
dipoleOriginActual = [0,0,zOffset]
dipoleOriginMirror = [0,0,-zOffset]
T = Theta/180*np.pi
P = Phi/180*np.pi
xComponent = np.sin(T) * np.cos(P)
yComponent = np.sin(T) * np.sin(P)
zComponent = np.cos(T)
normalizedPrimeActual = dipoleOriginActual[0] * xComponent + dipoleOriginActual[1] * yComponent + dipoleOriginActual[2] * zComponent
normalizedPrimeMirror = dipoleOriginMirror[0] * xComponent + dipoleOriginMirror[1] * yComponent + dipoleOriginMirror[2] * zComponent
phaseActual = np.exp(1j*k0*normalizedPrimeActual)
phaseMirror = -np.exp(1j*k0*normalizedPrimeMirror)
cos_psi = pol_vector[0]*xComponent + pol_vector[1]*yComponent + pol_vector[2]*zComponent
cos_psi = np.clip(cos_psi, -0.999, 0.999) # prevent divide by 0 error.
fieldIntensity = (np.cos(np.pi * ratio) - np.cos(np.pi * ratio * cos_psi)) / (cos_psi**2-1)
Ep_theta = pol_vector[0] * np.cos(T) * np.cos(P) + pol_vector[1] * np.cos(T) * np.sin(P) - pol_vector[2] * np.sin(T)
Ep_phi = -pol_vector[0] * np.sin(P) + pol_vector[1] * np.cos(P)
E_thetaActual = phaseActual * Ep_theta * fieldIntensity
E_phiActual = phaseActual * Ep_phi * fieldIntensity
E_thetaMirror = phaseMirror * Ep_theta * fieldIntensity
E_phiMirror = phaseMirror * Ep_phi * fieldIntensity
# Combine fields of actual dipole and mirror
E_theta = -1j * (E_thetaActual + E_thetaMirror)
E_phi = -1j * (E_phiActual + E_phiMirror )
if Theta > 90:
E_phi = 0
E_theta = 0
return [Phi, Theta, E_phi+1j*0, E_theta+1j*0] # return calculated farfield data point
[docs]
def createFarfield ( self, Type = "Dipole", intensity = 1, width = 80, antennaLength = 100, thetaMax = 180 ) :
"""
Creates a farfield of either a dipole or with a Gaussian distribution.
Parameters:
intensity (float): intensity modifier for Gaussian farfield.
width (float): maximum angle of the Gaussian farfield (in degrees).
antennaLength (float): length of the dipole section (oriented along the x-axis).
thetaMax (float): maximum value for theta (degree)
"""
self.intensity = intensity
self.width = width
self.antennaLength = antennaLength
self.printOut("creating farfield distribution with %s shape and wavelength %.2f" % (Type, self.Wavelength), self.log.FLOW_TYPE)
self.printOut("stepping: %.2f deg" % self.raster, self.log.DATA_TYPE, False)
if Type == self.__Type['GAUSS']:
self.printOut("1/e width: %.2f deg" % width, self.log.DATA_TYPE, False)
self.printOut("intensity: %.2f [au]" % intensity, self.log.DATA_TYPE, False)
elif Type == self.__Type['DIPOLE']:
self.printOut("Antenna-Length: %.2f mm" % antennaLength, self.log.DATA_TYPE, False)
self.FarField = [] # delete all data
for Phi in np.arange(0, 360 + self.raster, self.raster) : # loop all Phi (0 to 360. ==> 0deg is 360deg. two data points at 0/360
self.FarField.append([]) # add new phi-data list
f = self.FarField[-1]
for Theta in np.arange(0, thetaMax + self.raster, self.raster) : # loop all Theta (0 to 180. == no double point but 0 / 180 has df of 0)
if Type == self.__Type['GAUSS']:
farfieldDataPoint = self.calculateGaussianFarfieldPoint(Phi, Theta, intensity, width)
elif Type == self.__Type['DIPOLE']:
farfieldDataPoint = self.calculateDipoleFarfieldPoint(Phi, Theta, antennaLength)
else:
farfieldDataPoint = [Phi, Theta, 0+0j, 0+0j]
self.printOut("Farfield type '%s' doesn't exist. Farfield data point is set to 0." % Type, self.log.WARNING_TYPE)
f.append(farfieldDataPoint)
# Parameters which could be set in principle (so far not implemented)
# Coordinate system definitions (used to rotate and shift the data in 3D space) -- the following two parameters could be modified, that's why a de-rotation takes place below!!
self.OrientationX = [1, 0, 0]
self.OrientationZ = [0, 0, 1] # direction of propagation and direction with Theta = 0
self.Center = [0, 0, 0] # Center of the Far-field data in 3D space [mm]
self.PhaseCenter = [0, 0, 0] # Center of the phase-distribution, relative to the coordinate center [mm]
# Field data
self.Gain = 1 # Gain to boost the field if required
self.Phase = 0 # Phase of the file data [radian]
# Exitation data
self.Impedance = _Impedance # antenna impedance [Ohm]
# S-Parameters
self.S11 = 0 # S11
self.S21 = 1 # S21
# parameters resulting from the above settings / calculations
self.ReferenceDist = 1000 # default reference distance of the far-field pattern [mm] (1m set in df for integrating)
# calculate exitation current from surface integral over E^2 (Power = 0.5*I^2*Z = 0.5/Eta*Integral_over_|E|^2*df)
power = self.getPower() * 1e-6 # Reference distance is in mm not in m but power has units of W = kg m^2 / s^3
self.I0 = np.sqrt(2*power/self.Impedance)
self.Bandwidth = _Bandwidth # calculation benadwidth set to default value
# calculation of rotation matrices (nominal coordinate system is theta = 0 ==> positive z; Phi = 0 ==> positive x)
# calculate angles to rotate z-axis back to (0, 0, 1)
self.adjustFarfield()
self.Type = Type # set calculation type # Move this up to simplify code?
# update key-list if necessary
frequency = _C/self.Wavelength/1.e6
freqStr = "%.2f" % frequency
if freqStr not in self.freqDict:
self.freqDict[freqStr] = None
#self.plotFarField(figname = "generatedFarfield", show = True, logScale = False)
self.scaleI0()
#####################################
# combine far-field data
[docs]
def combineFarField ( self, FFData, dPhi = 5, dTheta = 5 ) :
"""
Combines far-field data of this instance with the once of an other instance and creates a new class instance.
Returns FarFieldDataClass with the resulting field.
Parameters:
FFData (FarFieldDataClass): with the far-field to combine. For phase variations use setPhase in before hand
dPhi (float): stepping in phi-direction of the new far-field [deg]
dTheta (float): stepping in theta-direction of the new far-field [deg]
Returns:
newFF (FarFieldDataClass): of the combined port
"""
self.printOut("Combining farfield data", self.log.FLOW_TYPE)
# check input data
# check wavelength
if self.Wavelength != FFData.Wavelength :
self.printOut("FarField data have different wavelength, can't proceed", self.log.ERROR_TYPE, False)
newFF = None
else:
self.printOut("setting up new farfield setup", self.log.DEBUG_TYPE, False)
# set up output far-field distribution
newFF = FarFieldData()
newFF.Wavelength = self.Wavelength
newFF.Impedance = self.Impedance
totalPower = (0.5*FFData.Impedance*FFData.I0**2) + (0.5*self.Impedance*self.I0**2)
newFF.I0 = np.sqrt(2*totalPower/newFF.Impedance)
newFF.FarField = []
# If fields come from same antenna arrangement in CST, far-field physically includes all S-parameters.
# Hence adding fields with phase-shifts directly results in the correct far-field distribution
phiRange = np.arange(0, 360 + dPhi, dPhi )
thetaRange = np.arange(0, 180 + dTheta, dTheta)
for phi in phiRange : # loop all new far-field points
self.printOut("calculating for phi = %.2f" % phi, self.log.DEBUG_TYPE, False, False)
dataLine = []
for theta in thetaRange :
data1 = self.getField(phi, theta)
data2 = FFData.getField(phi, theta)
dataLine.append([phi, theta, data1[2] + data2[2], data1[3] + data2[3]])
newFF.FarField.append(dataLine)
self.printOut("", self.log.DEBUG_TYPE, False)
power = newFF.getPower()
# implement new losses to the far-field data
newFF.S11 = (self.S11 + FFData.S11)/2
newFF.S21 = np.sqrt(power / (0.5*newFF.Impedance*newFF.I0**2))
if self.Type == FFData.Type :
newFF.Type = self.Type
else :
newFF.Type = self.__TYPE['NA']
return newFF
#####################################
# Methods to apply farField data
[docs]
def calculateCurrent ( self, ufo, add = True ) :
"""
Calculate current distribution induced by the farfield on a UFO file.
Parameters:
ufo (UFODataClass): If Type is _CURRENT the farfield inducted currents will be added
add (boolean): If the additional field and current distribution should be added to the existing
Returns:
ufoNew (UFODataClass): with the new data
"""
# create new ufo class and copy the header data
self.printOut("Calculating current on given UFO data class instance", self.log.FLOW_TYPE)
ufoNew = deepcopy(ufo)
ufoNew.Type = UFODataClass._CURRENT # new file type is _CURRENT
ufoNew.I0 = self.I0
ufoNew.Data = []
ufoNew.points = 0
if ufo.Type == UFODataClass._SURFACE : # don't add for surface data
add = False
# loop all data points within ufo and calculate the field at the surface points
self.printOut("calculating %d points on the surface given:" % len(ufo.Data), self.log.DEBUG_TYPE)
for idx in range(0, ufo.points) :
data = ufo.Data[idx]
E_, r_ , ph, th = self.getFieldPoint( data.v.x, data.v.y, data.v.z )
Er = [E_[0].real, E_[1].real, E_[2].real]
Ei = [E_[0].imag, E_[1].imag, E_[2].imag]
J_ = -2. * np.cross(np.asarray([data.n.x, data.n.y, data.n.z]), np.cross(r_, E_))/ _ETA # J = -2/eta * n_ x (r_ x E_)
Jr = [J_[0].real, J_[1].real, J_[2].real]
Ji = [J_[0].imag, J_[1].imag, J_[2].imag]
newD = UFODataClass.UFODataPoint()
newD.v = data.v
newD.n = data.n
if add :
newD.Er = UFODataClass._UFOVector(data.Er.x + Er[0], data.Er.y + Er[1], data.Er.z + Er[2])
newD.Ei = UFODataClass._UFOVector(data.Ei.x + Ei[0], data.Ei.y + Ei[1], data.Ei.z + Ei[2])
newD.Jr = UFODataClass._UFOVector(data.Jr.x + Jr[0], data.Jr.y + Jr[1], data.Jr.z + Jr[2])
newD.Ji = UFODataClass._UFOVector(data.Ji.x + Ji[0], data.Ji.y + Ji[1], data.Ji.z + Ji[2])
else :
newD.Er = UFODataClass._UFOVector( Er[0], Er[1], Er[2] )
newD.Ei = UFODataClass._UFOVector( Ei[0], Ei[1], Ei[2] )
newD.Jr = UFODataClass._UFOVector( Jr[0], Jr[1], Jr[2] )
newD.Ji = UFODataClass._UFOVector( Ji[0], Ji[1], Ji[2] )
E = np.array([newD.Er.x+1j*newD.Ei.x, newD.Er.y+1j*newD.Ei.y, newD.Er.z+1j*newD.Ei.z])
newD.I = np.dot(E, E.conj()).real
newD.df = data.df
ufoNew.Data.append(newD)
ufoNew.points = len(ufoNew.Data)
ufoNew.Data = (UFODataClass.UFODataPoint*ufoNew.points)(*ufoNew.Data) # establish c-type pointer array from data
self.printOut("Done", self.log.DEBUG_TYPE, False)
return ufoNew
#####################################
# Methods for plotting
[docs]
def plotFarField ( self, figname = "", show = False, logScale = True ) :
"""
Plotting the farfield distribution.
Parameters:
figname (string): filename of the created plot
show (boolean): flag if plot is shown (True) or safed to disk
"""
self.printOut("Plotting farfield distribution", self.log.FLOW_TYPE)
phiRange = range(0, 360+self.raster, self.raster)
thetaRange = range(0, 180+self.raster, self.raster)
intens = np.zeros((len(phiRange), len(thetaRange)))
for p_idx, phi in enumerate(phiRange) :
for t_idx, theta in enumerate(thetaRange) :
f = self.getField(phi, theta) #phi, th, Ep, Et
if logScale:
intens[p_idx, t_idx] = 10*np.log10(np.absolute(f[2])**2 + np.absolute(f[3])**2 + .0000001)+60
else:
intens[p_idx, t_idx] = np.sqrt(np.absolute(f[2])**2 + np.absolute(f[3])**2 )
Theta = np.array(thetaRange)/180*np.pi
Phi = np.array(phiRange)/180*np.pi
THETA, PHI = np.meshgrid(Theta, Phi)
X = intens * np.sin(THETA) * np.cos(PHI)
Y = intens * np.sin(THETA) * np.sin(PHI)
Z = intens * np.cos(THETA)
fig = plt.figure()
ax = plt.axes(projection='3d')
ax.xaxis.pane.fill = False
ax.xaxis.pane.set_edgecolor('white')
ax.yaxis.pane.fill = False
ax.yaxis.pane.set_edgecolor('white')
ax.zaxis.pane.fill = False
ax.zaxis.pane.set_edgecolor('white')
ax.grid(False)
ax.set_zlim(-np.max(Z), np.max(Z))
ax.set_xlim(-np.max(Z), np.max(Z))
ax.set_ylim(-np.max(Z), np.max(Z))
ax.set_aspect('equal')
# plt.axis('off')
plot = ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=plt.get_cmap('jet'), linewidth=0, antialiased=False, alpha=0.5)
if figname != "" :
fig.savefig(figname)
if show :
plt.show()
else :
plt.show()
plt.close()
[docs]
def plotFarFieldCut ( self, phi = 0, polAngle = 0, dB = True, figname = "", show = False ) :
"""
Plotting the farfield distribution as cut on a given phi-angle.
Parameters:
phi (float or list of floats): phi angle of the cut
polAngle (float): Phi-angle of the polarization to look at
dB (boolean): Indicating if a dB scale (True) or a lineare scale is used
figname (string): Filename of the created plot
show (boolean): Flag if plot is shown (True) or safed to disk
"""
self.printOut("Plotting farfield distribution as cut", self.log.FLOW_TYPE)
if isinstance(phi, float) or isinstance(phi, int) :
phi = [phi]
colors = getColors(len(phi))
# set figure details and labels
fig, ax = plt.subplots()
# extract data
maxT = self.FarField[0][-1][1]
thetaRange = np.arange(-maxT, maxT+self.raster, self.raster)
intens = []
for p in phi :
intens.append(np.zeros ( len(thetaRange) ) )
iTotal = np.zeros ( len(thetaRange) )
diffPhi = np.pi*(polAngle)/180
for idx, theta in enumerate(thetaRange) :
if theta < 0 :
f = self.getField(p+180, -theta)
else :
f = self.getField(p, theta)
f = np.asarray(f)
intens[-1][idx] = np.absolute( f[2] * np.sin(diffPhi) + f[3] * np.cos(diffPhi) )**2
iTotal[idx] = np.absolute( f[2]*f[2].conj() + f[3]*f[3].conj() )
if dB :
intens[-1] = intens[-1]/np.amax(iTotal)
intens[-1][intens[-1] < 1e-99] = 1e-99
intens[-1] = np.log10(intens[-1])*10
if dB :
ax.set_ylabel("Intensity [dB]")
else :
ax.set_ylabel("Intensity [au]")
unit = "degree"
if maxT < self.raster :
thetaRange *= 60
unit = "minutes of arc"
if maxT * 60 < self.raster :
thetaRange *= 60
unit = "seconds of arc"
ax.set_xlabel("Theta angle [%s]" % unit)
ax.set_title ("%.2fGHz : Field intensity cuts)" % (_C/self.Wavelength/1e6))
for idx, p in enumerate(phi) :
ax.plot(thetaRange[intens[idx] > -100], intens[idx][intens[idx] > -100], label = "phi = %.1f°" % p, color = colors[idx], linestyle='solid')
fig.legend(loc='upper left', bbox_to_anchor=(0.65, 0.85))
# check safing the figure
if figname != "" :
fig.savefig(figname)
if show :
plt.show()
else :
plt.show()
plt.close()
[docs]
def plotFarField2d ( self, figname = "", show = False, rotate = False, upside = False, vmin = -30 ) :
"""
Plotting the farfield distribution as 2d-false color plot including rotation.
Parameters:
figname (string): Filename of the created plot.
show (boolean): Flag if plot is shown (True) or saved to disk.
rotate (boolean): Flag indicating if rotation angels are included or not
upside (boolean): Flip z-axis (has only effect if rotate is set to True)
vmin (float): Indicating lower bound for dB scale
"""
self.printOut("Plotting farfield distribution in 2d", self.log.FLOW_TYPE)
# set figure details and labels
fig, ax = plt.subplots()
ax.set_frame_on(False)
ax.axis('off')
ax.set_aspect('equal')
ax.set_xlabel("Radial: Theta; axial : Phi")
ax.set_title ("%.2fGHz : Field intensity [dB]\n(contours are at -3, -10, and -20dB)" % (_C/self.Wavelength/1e6))
# extract data
maxT = self.FarField[0][-1][1]
thetaRange = np.arange(-maxT, maxT+self.raster, self.raster)
intens = np.zeros((len(thetaRange) **2, 3))
idx = 0
for x in thetaRange :
for y in thetaRange :
phi = np.arctan2(y, x)*180/np.pi
theta = np.sqrt(x**2 + y**2)
if theta < maxT :
f = self.getField(phi, theta)
else :
f = [0,0,0,0]
intens[idx, :] = np.asarray([x, y, np.absolute(f[2])**2. + np.absolute(f[3])**2])
idx += 1
# set x-y plane and scale z
x = intens[:, 0]
y = intens[:, 1]
z = intens[:, 2]
z = z / max(z)
x = x[z != 0]
y = y[z != 0]
z = z[z != 0]
z = 10.*np.log10(z)
# re-grid x-y-z
xi, yi = np.mgrid[min(x):max(x):_GridSize*1j, min(y):max(y):_GridSize*1j]
zi = griddata((x, y), z, (xi, yi), method = 'linear')
# plot data
cf = ax.pcolormesh(xi, yi, zi, vmin= vmin, vmax= 0., cmap = plt.get_cmap('rainbow'))
ax.contour(xi, yi, zi, levels=[-20, -10, -3], linewidths = 1.5, colors = 'w')
fig.colorbar(cf)
# plot phi-scale with labels
# setup phi-scale circle
circleX = []
circleY = []
for p in range(0, 340) :
xp = np.cos(p/180.*np.pi)*maxT
yp = np.sin(p/180.*np.pi)*maxT
circleX.append( xp )
circleY.append( yp )
if p == int(p/10)*10 :
circleX.append(xp*0.97)
circleY.append(yp*0.97)
circleX.append(xp*1.03)
circleY.append(yp*1.03)
circleX.append(xp)
circleY.append(yp)
circleX = np.asarray(circleX)
circleY = np.asarray(circleY)
ax.plot(circleX, circleY, 'k')
ax.text(0., maxT*1.05, "phi= 90°", ha = 'center', va = 'bottom')
ax.text(0., -maxT*1.05, "phi= -90°", ha = 'center', va = 'top' )
ax.text(-maxT*1.05, 0., "phi= 180°", ha = 'right', va = 'center', rotation = 90)
ax.text( maxT*1.05, 0., "phi= 0°", ha = 'left', va = 'center', rotation = -90)
# plot theta = 90deg circle to indicate forward hemisphere
if maxT > 90. :
circleX = circleX/maxT*90.
circleY = circleY/maxT*90.
ax.plot(circleX, circleY, linewidth = 0.8)
# plot theta scale with labels
# setup theta scale
lineX = []
lineY = []
for t in range(0, int(maxT), 10) :
lineX.append(t )
lineY.append(0.)
if int(t/90)*90 == t :
lineX.append(t )
lineY.append( 0.03*maxT)
lineX.append(t )
lineY.append(-0.03*maxT)
else:
lineX.append(t )
lineY.append( 0.01*maxT)
lineX.append(t )
lineY.append(-0.01*maxT)
lineX.append(t )
lineY.append(0.)
ax.plot(lineX, lineY, 'k')
ax.text( maxT, -0.05*maxT, "theta = %d°" % int(maxT), ha = 'center', va = 'top', rotation = -90)
ax.text( 90., -0.05*maxT, "theta = 90°" , ha = 'center', va = 'top', rotation = -90)
# check safing the figure
if figname != "" :
fig.savefig(figname)
if show :
plt.show()
else :
plt.show()
plt.close()
[docs]
def makeMovie ( self, movName = "farfield.mov", folder = "" ) :
"""
Creates a movie on all far-field data in the folder set (if a folder was set).
Parameters:
movName (string): name of the movie
folder (string): folder to store the individual plots and the movie
Returns:
boolReturn (boolean): flag to indicate whether a movie was created.
"""
self.printOut("Creating a movie on all far-field data of the folder set", self.log.FLOW_TYPE)
if len(self.freqDict) < 2 :
self.printOut("not enough far-field data or folder not set", self.log.ERROR_TYPE, False)
boolReturn = False
else:
for idx, frequency in enumerate(sorted(self.freqDict)) :
self.setFrequency(float(frequency))
self.plotFarField2d("%sFarField_%dMHz.png" % (folder, int(float(frequency)*1000+0.5)), False)
makeMovie(folder + "FarField_*MHz.png", folder + movName, len(self.freqDict)/2. )
boolReturn = True
return boolReturn
[docs]
def plotSParam ( self, figname = "", show = False ) :
"""
Plots the s-parameter data of the folder as set.
Parameters:
figname (string): Filename of the created plot.
show (boolean): Flag if plot is shown (True) or saved to disk.
"""
# initialize s-parameter data class for plot
self.printOut("Plotting the efficiency and if available S-parameter data of the farfields in the folder", self.log.FLOW_TYPE)
if (len(self.freqDict) < 2) :
self.printOut("not enough farfield data or folder not set", self.log.ERROR_TYPE, False)
freq = np.zeros(len(self.freqDict))
s11 = np.zeros(len(self.freqDict))
rad = np.zeros(len(self.freqDict))
for idx, frequency in enumerate(sorted(self.freqDict)) :
self.setFrequency(float(frequency))
freq[idx] = float(frequency)
s11 [idx] = np.log10(np.absolute(self.S11))*20.
rad [idx] = np.log10(np.absolute(self.S21))*20.
ledgends = []
fig, ax = plt.subplots()
for label in ax.get_yticklabels():
label.set_color("blue")
ax2 = ax.twinx()
for label in ax2.get_yticklabels():
label.set_color("red")
data = self.cst.getItemData("S1,1", np.amin(freq), np.amax(freq))
if data == None :
ax.plot(freq, s11, "b-")
ledgends.append("Returned power from far-field")
else :
ax.plot(freq, s11, "bo")
ledgends.append("Returned power from far-field")
ax.plot(data[0], np.log10(np.absolute(data[1]))*20., "b-")
ledgends.append("$|S11|^2$ from CST")
data = self.cst.getItemData("Tot. Efficiency", np.amin(freq), np.amax(freq))
if data == None :
ax2.plot(freq, rad, "r-")
ledgends.append("total efficiency from far-field")
else:
ax2.plot(freq, rad, "r+")
ledgends.append("total efficiency from far-field")
ax2.plot(data[0], np.log10(np.absolute(data[1]))*10., "r-")
ledgends.append("total efficiency from CST")
ax2.set_ylabel("Efficiency [dB]", color="red")
ax.set_ylabel("$|S11|^2$ / returned power [dB]", color = "blue")
ax.set_xlabel("Frequency [GHz]")
fig.legend(ledgends, loc='lower left', bbox_to_anchor=(0.1, 0.1))
if figname != "" :
fig.savefig(figname)
if show :
plt.show()
else :
plt.show()
plt.close()
[docs]
def plotOpeningAngle ( self, figname = "", show = False ) :
"""
Plots the average opening angle over frequency (if a folder is set).
Parameters:
figname (string): Filename of the created plot.
show (boolean): Flag if plot is shown (True) or saved to disk.
"""
self.printOut("Plotting the opening angle of the farfields in the folder", self.log.FLOW_TYPE)
if (len(self.freqDict) < 2) :
self.printOut("not enough farfield data or folder not set", self.log.ERROR_TYPE, False)
# extract 3, 10 and 20 dB opening angle
openingAngle = []
freq = []
for idx, frequency in enumerate(sorted(self.freqDict)) :
self.setFrequency(float(frequency))
freq.append(float(frequency))
av3 = []
av10 = []
av20 = []
for phi in range (0, 360, 1) :
a3 = False
a10 = False
a20 = False
for th in range (0, 180, 1) :
p, t, ep, et = self.getField(phi, th)
I = np.absolute(ep)**2. + np.absolute(et)**2.
if th == 0:
ref = I
old = I/ref
v = I/ref
if (v < 0.5) and (not a3) :
av3.append(th + (old-0.5)/(old-I) - 1.)
a3 = True
if (v < 0.1) and (not a10) :
av10.append(th + (old-0.1)/(old-I) - 1.)
a10 = True
if (v < 0.01) and (not a20) :
av20.append(th + (old-0.01)/(old-I) - 1.)
a20 = True
old = v
if len(av20) == 0:
av20.append(180.)
if len(av10) == 0:
av10.append(180.)
if len(av3) == 0:
av3.append(180.)
openingAngle.append([
np.amin(np.asarray(av3)),
np.amax(np.asarray(av3)),
np.average(np.asarray(av3)),
np.amin(np.asarray(av10)),
np.amax(np.asarray(av10)),
np.average(np.asarray(av10)),
np.amin(np.asarray(av20)),
np.amax(np.asarray(av20)),
np.average(np.asarray(av20)),
] )
openingAngle = np.asarray(openingAngle)
fig, ax = plt.subplots()
ax.plot(freq, openingAngle[:, 0], 'k--')
ax.plot(freq, openingAngle[:, 1], 'k:' )
ax.plot(freq, openingAngle[:, 2], 'k-' )
ax.plot(freq, openingAngle[:, 3], 'r--')
ax.plot(freq, openingAngle[:, 4], 'r:' )
ax.plot(freq, openingAngle[:, 5], 'r-' )
ax.plot(freq, openingAngle[:, 6], 'b--')
ax.plot(freq, openingAngle[:, 7], 'b:' )
ax.plot(freq, openingAngle[:, 8], 'b-' )
fig.legend(['-3dB (min)', '-3dB (max)', '-3dB (aver)',
'-10dB (min)', '-10dB (max)', '-10dB (aver)',
'-20dB (min)', '-20dB (max)', '-20dB (aver)' ], loc='upper left', bbox_to_anchor=(0.75, 0.85))
# establish labling
ax.set_xlabel("Frequency [GHz]")
ax.set_ylabel("opening angle [deg] ( half-angle)")
ax.set_title ("opening angle over frequency for several power-level" )
# safe of show figure?
if figname != "" :
fig.savefig(figname)
if show :
plt.show()
else :
plt.show()
plt.close()
#####################################
# Methods to change coordinate system
[docs]
def movePhaseCenter ( self, x, y, z ) :
"""
Moves phase center (rotation-point of the far-field distribution) to point (x, y, z)
Parameters:
x (float): x coordinate of new point
y (float): y coordinate of new point
z (float): z coordinate of new point
"""
self.printOut("Moving phase center to (%.2f, %.2f, %.2f)" %(x, y, z), self.log.DEBUG_TYPE)
self.PhaseCenter = [x, y, z]
[docs]
def moveCenter ( self, x, y, z ) :
"""
Moves the FarField coordinate system to point (x, y, z)
Parameters:
x (float): x coordinate of new point
y (float): y coordinate of new point
z (float): z coordinate of new point
"""
self.printOut("Moving center to (%.2f, %.2f, %.2f)" %(x, y, z), self.log.DEBUG_TYPE)
self.Center = [ x, y, z]
[docs]
def resetRotation ( self ) :
"""
Sets coordinate system rotation to identity matrix (no rotation)
"""
self.printOut("Reseting rotation data", self.log.DEBUG_TYPE)
# Modify coordinate vectors with stored rotation data
self.OrientationX = np.dot( self.invRot, self.OrientationX )
self.OrientationZ = np.dot( self.invRot, self.OrientationZ )
self.Rot = np.identity(3) # Rotation matrix for the output data
self.invRot = np.linalg.inv(self.Rot) # inverse rotation
[docs]
def rotate ( self, alpha_x, alpha_y, alpha_z ) :
"""
Rotates the actual coordinate system relative to x-, y- and z-axis
Parameters:
alpha_x (float): angle of rotation relative to x axis [Deg]
alpha_y (float): angle of rotation relative to y axis [Deg]
alpha_z (float): angle of rotation relative to z axis [Deg]
"""
self.printOut("Rotating coordinates by alpha_x = %.2f, alpha_y = %.2f, alpha_z = %.2f" %(alpha_x, alpha_y, alpha_z), self.log.DEBUG_TYPE)
alpha_x = alpha_x/180.*np.pi
alpha_y = alpha_y/180.*np.pi
alpha_z = alpha_z/180.*np.pi
# set up rotation matrices
Rotx = np.array( [ [ 1., 0., 0.],
[ 0., np.cos(alpha_x), -np.sin(alpha_x)],
[ 0., np.sin(alpha_x), np.cos(alpha_x)] ])
Roty = np.array( [ [ np.cos(alpha_y), 0., np.sin(alpha_y)],
[ 0., 1., 0.],
[ -np.sin(alpha_y), 0., np.cos(alpha_y)] ])
Rotz = np.array( [ [ np.cos(alpha_z), -np.sin(alpha_z), 0.],
[ np.sin(alpha_z), np.cos(alpha_z), 0.],
[ 0., 0., 1.] ])
# and combine them with each other and
Rot_n = np.dot(Rotz, np.dot(Roty, Rotx))
# with the existing rotation matrix
self.Rot = np.dot( Rot_n, self.Rot)
# determin the inverse rotation
self.invRot = np.linalg.inv(self.Rot)
# Modify coordinate vectors with stored rotation data
self.OrientationX = np.dot( Rot_n, self.OrientationX )
self.OrientationZ = np.dot( Rot_n, self.OrientationZ )
#####################################
# Methods to work with the data
[docs]
def getPower ( self ) :
"""
Calculates the actual power in the far-field distribution and compares it with the excitation power via the excitation current.
Returns:
p (float): Far-Field power [W]
"""
stepP = np.pi/180
stepT = np.pi/180
p = 0
for Phi in range(360) :
for Theta in range(180) :
data = self.getField(Phi, Theta)
p += (np.absolute(data[2])**2 + np.absolute(data[3])**2) * np.sin( data[1] * np.pi / 180 ) * stepP * stepT * self.ReferenceDist**2.
p = p/_ETA/2
self.printOut( "Power = %.2f W, relative power = %.2f %%" %( p, 100 * p / ( 0.5 * self.Impedance * self.I0**2 ) ), self.log.DATA_TYPE )
return p
[docs]
def scaleI0 ( self, newI0 = _I0iso ) :
"""
Scales the field-amplitudes such that they fullfill the new excitation current without changing the input impedance.
Parameters:
newI0 (float): new excitation current
"""
self.printOut("Scaling I0 from %.2f A to %.2f A" % (self.I0, newI0), self.log.DEBUG_TYPE)
for Phi in self.FarField :
for data in Phi :
data[2] = newI0*data[2]/self.I0*self.Gain
data[3] = newI0*data[3]/self.I0*self.Gain
self.Gain = 1.
self.I0 = newI0
[docs]
def setGain ( self, gain, relative = True ) : # set Gain parameter
"""
Sets the gain of the system, which means the interpolated field is multiplied with the scalar factor gain.
Parameters:
gain (float): Factor to multiply the field with (set to zero gain, means correct back)
relative (boolean): Correct the I0 value (True) or introduce a loss (False)
"""
self.printOut("Setting gain to %.2f" % gain, self.log.DEBUG_TYPE)
if relative :
oldGain = self.Gain
self.Gain *= gain
if self.Gain == 0. :
self.Gain = 1.
self.I0 /= oldGain
else :
self.I0 *= gain
else :
self.S21 *= np.sqrt(gain)
for fp in self.FarField :
for f in fp :
f[2] *= gain
f[3] *= gain
[docs]
def setPhase ( self, phase, absolute = False ) :
"""
Sets the phase of the system, which means that the field points are shifted in phase by this value.
Parameters:
phase (float): new phase [deg]
absolute (boolean): phase is added to the actual value (False) or taken as new absolute phase value (True)
"""
self.printOut("Setting phase to %.2f radians" % phase, self.log.DEBUG_TYPE)
if absolute : self.Phase = 0.
self.Phase += phase
[docs]
def getField ( self, phi, th, interpolate = False ) :
"""
Method to interpolate a point within the loaded FarField data array.
Parameters:
phi (float): Angle Phi for which the data is interpolated
th (float): Angle Theta for which the data is interpolated
interpolate (bool): Whether or not the point should be interpolated or nearest
Returns:
farfieldPoint (list): The point on the farfield at phi,th in format [phi, th, Ephi, Eth]
"""
def bilinear_interpolation(phi, th, surrPoints, idx):
"""
Perform bilinear interpolation of point phi, th using surrounding surrPoints
Parameters:
phi (flaot): Phi coordinate of interpolation point
th (flaot): Theta coordinate of interpolation point
surrPoints (list): List of coordinates for the four points surrounding interpolation point
idx (int): Index of the value (2 for phi component, 3 for theta component)
Returns:
result (float): Interpolated value idx (phi or theta component) at phi,th
"""
phiStart, thStart = surrPoints[0][:2]
phiEnd, thEnd = surrPoints[3][:2]
result = (
surrPoints[0][idx] * (phiEnd - phi ) * (thEnd - th ) +
surrPoints[1][idx] * (phi - phiStart) * (thEnd - th ) +
surrPoints[2][idx] * (phiEnd - phi ) * (th - thStart) +
surrPoints[3][idx] * (phi - phiStart) * (th - thStart)
) / ((phiEnd - phiStart) * (thEnd - thStart))
return result
phi = phi % 360 # 360deg = 0deg
th = th % 180 # 180deg = 0deg
# Skipping interpolation is still experimental.
#if self.Type == "DIPOLE":
# farfieldPoint = self.calculateDipoleFarfieldPoint(phi, th, self.antennaLength)
#elif self.Type == "GAUSS":
# farfieldPoint = self.calculateGaussianFarfieldPoint(phi, th, self.intensity, self.width)
if interpolate:
phStart = int(phi / self.raster) * self.raster # phi coordinate before interpolation point
thStart = int(th / self.raster) * self.raster # theta coordinate before interpolation point
surrPoints = [] # list of the four points surrounding phi,th
surrPoints.append(self.getField(phStart , thStart )) # add the four points to the list
surrPoints.append(self.getField(phStart + self.raster, thStart ))
surrPoints.append(self.getField(phStart , thStart + self.raster))
surrPoints.append(self.getField(phStart + self.raster, thStart + self.raster))
phiInterpolate = bilinear_interpolation(phi, th, surrPoints, 2) # perform interpolation for index 2 (phi component of farfield)
thetaInterpolate = bilinear_interpolation(phi, th, surrPoints, 3) # perform interpolation for index 3 (theta component of farfield)
farfieldPoint = [phi, th, phiInterpolate, thetaInterpolate]
else:
phIdx = round(phi / self.raster) # array index in phi-direction
thIdx = round(th / self.raster) # array index in theta-direction
Ep = self.Gain*self.FarField[phIdx][thIdx][2]*np.exp(1j*self.Phase/180.*np.pi) # multiply with complex gain (Gain and Phase) and the normalization
Et = self.Gain*self.FarField[phIdx][thIdx][3]*np.exp(1j*self.Phase/180.*np.pi)
farfieldPoint = [phi, th, Ep, Et]
return farfieldPoint
[docs]
def getFieldPoint ( self, xp, yp, zp ) :
"""
Calculate E-field at position (xp, yp, zp).
Parameters:
xp (float): x coordinate
yp (float): y coordinate
zp (float): z coordinate
Returns:
complex 3d E-vector (?)
normalized pointing vector (from FarField center-coordinate to surface point)
"""
# calculate position within the FarField frame:
pos_in = np.array([xp, yp, zp]) # set desired position in far-field coordinate system
pos_in = pos_in - self.Center # remove center position
pos = np.dot(self.invRot, pos_in) # inverse rotate coordinate system
pos = pos - self.PhaseCenter # substract phase center
# calculate spherical coordinates
r = np.sqrt(pos[0]**2. + pos[1]**2. + pos[2]**2.) # distance between phase center and surface point
thRad = np.arccos(pos[2]/r) # calculate Theta in radians
th = thRad * 180. / np.pi # and in degree
phRad = np.arctan2(pos[1], pos[0]) # calculate Phi in radians
ph = phRad * 180. / np.pi # and in degree
ph = ph % 360 # adjust angle to 0 to 360.deg
# determin field at coordinate positions
#f = self.getField(ph, th) # interpolate field position
f = self.getField(ph, th, interpolate = True) # interpolate field position
# adjust phase and amplitude to a distance r taking the reference distance into account
Ephi = self.ReferenceDist * f[2] /r * np.exp(-1j*np.pi*2.*(r-self.ReferenceDist)/self.Wavelength)
Eth = self.ReferenceDist * f[3] /r * np.exp(-1j*np.pi*2.*(r-self.ReferenceDist)/self.Wavelength)
# transform field to cartesian coordinates
E = np.array( [ Eth * np.cos( thRad ) * np.cos( phRad ) - Ephi * np.sin( phRad ),
Eth * np.cos( thRad ) * np.sin( phRad ) + Ephi * np.cos( phRad ),
-Eth * np.sin( thRad )
] )
E = np.dot(self.Rot, E)
return E, pos_in/r, ph, th
[docs]
def writeSphere ( self, name, radius = 100., thetaR = [0, 90] ) :
"""
Writes a cartesian coordinate sphere to an ascii data file.
Parameters:
name (string): output file name
radius (float): radius of the sphere [mm]
"""
self.printOut("Writing sphere data to file %s" % name, self.log.FLOW_TYPE)
self.printOut("using radius %.2f mm" % radius, self.log.DATA_TYPE, False)
self.printOut("using Teta-Range [%.2f, %.2f] deg" % thetaR, self.log.DATA_TYPE, False)
fh = open(name, 'w') # open file for writing
for th in range(thetaR[0], thetaR[1]) :
for ph in range(360):
if (th == 0) and (ph > 0) : continue
thRad = th/180.*np.pi
phRad = ph/180.*np.pi
x = radius*np.sin(thRad)*np.cos(phRad)
y = radius*np.sin(thRad)*np.sin(phRad)
z = radius*np.cos(thRad)
E, r, phn,thn = self.getFieldPoint(x, y, z)
I = np.absolute(E[0])**2. + np.absolute(E[1])**2. + np.absolute(E[2])**2.
fh.write("%.4f %.4f %.4f %.4f %.4f %.4f %.4f %.6f %.6f %.6f %.6f %.6f %.6f %.6f %.4f\n"
% (x, y, z, th, ph, thn, phn, E[0].real, E[1].real, E[2].real, E[0].imag, E[1].imag, E[2].imag, I, radius**2*np.sin(thRad)))
[docs]
def estimatePhaseCenter ( self ) :
"""
Try to estimate the phase center of the given far field distribution
The result is stored to PhaseCenter
"""
def getRMS(xp, yp, zp) :
radius = 101.234235 * self.Wavelength # 101.2345 chosen to avoid gridding issues.
self.moveCenter(xp, yp, zp)
SphX = 0.
SphXX = 0.
SphY = 0.
SphYY = 0.
n = 0
for th in range (1, 92, self.raster) :
thRad = th/180.*np.pi
for ph in range(0, 360, 20) :
phRad = ph/180.*np.pi
x = radius*np.sin(thRad)*np.cos(phRad)
y = radius*np.sin(thRad)*np.sin(phRad)
z = radius*np.cos(thRad)
E, r, phn,thn = self.getFieldPoint(x, y, z)
phaseX = np.arctan2(E[0].imag, E[0].real)
phaseY = np.arctan2(E[1].imag, E[1].real)
SphX += phaseX
SphXX += phaseX*phaseX
SphY += phaseY
SphYY += phaseY*phaseY
n += 1
mx = SphX/n
my = SphY/n
rmsX = ( SphXX - 2*mx*SphX + n*mx*mx )/n
rmsY = ( SphYY - 2*my*SphY + n*my*my )/n
return np.sqrt(rmsX+rmsY)
self.printOut("Estimating phase-center of CST-file -- please be patient .... ", self.log.FLOW_TYPE)
minRMS = 1.E10
minX = 0
minY = 0
minZ = 0
xp = 0
yp = 0
zp = 0
for i in range(2, -1, -1) :
step = 10.**i
xc = minX
yc = minY
zc = minZ
for xpp in range(-6, 7, 1) :
xp = xpp*step + xc
for ypp in range(-6, 7, 1) :
yp = ypp*step + yc
for zpp in range (-6, 7, 1) :
zp = zpp*step + zc
rms = getRMS(xp, yp, zp)
if rms < minRMS :
minRMS = rms
minX = xp
minY = yp
minZ = zp
self.printOut("with step %.2f ==> phase-center at (%.2f, %.2f, %.2f); RMS = %.2f" %(step, minX, minY, minZ, minRMS), self.log.DEBUG_TYPE, False, False)
self.printOut("", self.log.DEBUG_TYPE, False)
self.PhaseCenter = [minX, minY, minZ]
self.Center = [ 0., 0., 0.]
#############################################################################
# Main function; -- for test-purposes
#
#############################################################################
if __name__=="__main__":
ff = FarFieldData(100)
#ff.setFolder(root_path_string+"/Input_Files/CST_Files/", wildcard = "?.??_GHz.ffs")
ff.createGaussianFarField()
#ff.createDipoleFarField()
ff.plotFarField(show = True)#, logScale = False)
#ff.plotFarField2d(show = True)