#!/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
#############################################################################
#
# Provide calculation tools on Gaussian beam propagation in vacuum (refractive index n == 1 ).
# In principle all input data must only have the same length unit and the output data in consequence have the same unit. Hence coice of the length unit is in principle free.
# The UFO modules use per convention millimeters as length unit. And there are functions and classes expecting input millimeters. Hence using other length units as millimeters
# is possible for this class, but may requires unit conversion when using other classes or functions of the package.
#
#
##############################################################################
# - clean-up code (frequently and regularly, please)
# - improve readability and add more comments (frequently and regularly)
#
# single task-changes / add-ons
#
#
#############################################################################
# import libs
#############################################################################
import numpy as np
from copy import copy, deepcopy
import LogClass
#############################################################################
# constants
#############################################################################
_CLASSNAME = "UFO-Gauss-Calc"
log = LogClass.LogClass(_CLASSNAME) # Logging class used for all in and output messages
#############################################################################
# Methods
#############################################################################
[docs]
def calcDist2PhaseCenter ( w, d, wavelength) :
"""
calculates the distance to the phase-center from a given waist w and an distance d from the waist.
Parameters:
w (float): waist [mm]
d (float): distance from the waist w [mm]
wavelength (float): wavelength for the calculation [mm]
Returns:
float, distance r to the nominal phase-center [mm]. The interpretation is as follows: r denotes the radius of the phase-front in a distance d from the waist w.
or None for an input distance d == 0
"""
zr = np.pi * w * w / wavelength # calculate the Rayleigh range
if d == 0. : return None # check for a distance of zero ==> return None
return d + zr*zr / d # return radius of curvature
[docs]
def calcConicData ( w1, w2, d1, d2, wavelength, reflectionAngle, verbose = True, log = log ) :
"""
Calculate data of a ellipsoide or hyperbola for the Gaussian beam parameters given. A distance d1 or d2 of zero is interpreted as 10^-6 mm distance to avoid numerical problems
==> Parabolas get very large ellipses.
Parameters:
w1 (float): waist-size 1 [mm]
w2 (float): waist-size 2 [mm]
d1 (float): distance of waist w1 from the conical surface [mm]
d2 (float): distance of waist w2 from the conical surface [mm]
wavelength (float): wavelength used for the calculation [mm]
reflectionAngle (float): reflection angle in degree
verbose (bool): does nothing
Returns:
2d conic section data in x-y plane. Origin of the coordinate system is the center of the conical form. The function returns all data including redundent information
a (float): semi-major axis [mm]
b (float): semi-minor axis (negative for hyperbola) [mm]
e (float): excentricity
c (float): linear eccentrisity [mm]
p (float): size factor [mm]
ref.point-x (float): reflection point x-coordinate at the conical section
ref.point-y (float): reflection point y-coordinate at the conical section
tangentangle (float): tangent angle at the reflection point (angle between major axis and tangent [degree]
focal-length (float): focal length of the conic section at the reflection point [mm]
ref.angle (float): reflection angle [degree]
r1 (float): distance of waitst w1 to its phase center [mm]
r2 (float): distance of waitst w2 to its phase center [mm]
"""
log.printOut("Calculating conic data for :", log.DATA_TYPE)
log.printOut(" w1 = %.2f mm" % ( w1 ), log.DATA_TYPE, False)
log.printOut(" d1 = %.2f mm" % ( d1 ), log.DATA_TYPE, False)
log.printOut(" w2 = %.2f mm" % ( w2 ), log.DATA_TYPE, False)
log.printOut(" d2 = %.2f mm" % ( d2 ), log.DATA_TYPE, False)
log.printOut(" ref. Angle = %.2f deg" % ( reflectionAngle ), log.DATA_TYPE, False)
log.printOut(" wavelength = %.2f mm" % ( wavelength ), log.DATA_TYPE, False)
# calculare Rayleigh ranges of the beams
zr1 = np.pi * w1 * w1 / wavelength
zr2 = np.pi * w2 * w2 / wavelength
# check for zero distance to waist and exchange with low number to avoid numerical problems (could also be a seperate code for parabolas)
if d1 == 0. : d1 = 1e-3
if d2 == 0. : d2 = 1e-3
# calculate distances to phase center
r1 = d1 + zr1*zr1 / d1
r2 = d2 + zr2*zr2 / d2
# focal parameter p
p = r1 * r2 * ( 1 + np.cos( np.pi * reflectionAngle / 180 ) ) / ( r1 + r2 )
# eccentricity e
e = np.sqrt( 1 + ( r2 / ( r1+ r2 ) - 1 ) * 2 * p / r1 )
# semi-axis a and b
bb = p * p / ( 1 - e * e ) # bb < 0 for a hyperbola!!
aa = bb * bb / ( p * p )
# linear eccentricity c
c = np.sqrt( aa - bb )
# get reflection point in internal coordinates (x-z plane)
# avoid negative values in sqare-root (angle = 90°) ... only parabola with r1 or r2 = inf
if ( r1 * r1 - ( p - r1 ) * ( p - r1 ) / ( e * e ) ) < 0. :
center = [ 0., 0., ( p - r1 ) / e ]
else :
center = [ np.sqrt( r1 * r1 - ( p - r1 ) * ( p - r1 ) / ( e * e ) ), 0., ( p - r1 ) / e ]
turn = False # flip structure (mirror on y-axis) for hyperbola on negative x-values
tangens = np.arctan2( ( aa * center[0] ), - r1 / np.abs(r1) * bb * ( center[2] + c ) )
b = np.sqrt(np.abs(bb))
if bb < 0 : b = b*1j
log.printOut(" Results:", log.DATA_TYPE, False)
log.printOut(" a = %.2f mm" % ( np.sqrt(aa) ), log.DATA_TYPE, False)
if b.imag != 0. :
log.printOut(" b = %.2f mm (imaginary)" % ( b.imag ), log.DATA_TYPE, False)
else :
log.printOut(" b = %.2f mm" % ( b ), log.DATA_TYPE, False)
log.printOut(" e = %.3f" % ( e ), log.DATA_TYPE, False)
log.printOut(" c = %.3f mm" % ( c ), log.DATA_TYPE, False)
log.printOut(" p = %.3f mm" % ( p ), log.DATA_TYPE, False)
log.printOut(" center = ( %.2f mm, %.2f mm )" % ( center[0], center[2] ), log.DATA_TYPE, False)
log.printOut(" angle = %.2f deg" % ( tangens*180./np.pi ), log.DATA_TYPE, False)
log.printOut(" r1 = %.2f mm" % ( r1 ), log.DATA_TYPE, False)
log.printOut(" r2 = %.2f mm" % ( r2 ), log.DATA_TYPE, False)
return np.sqrt(aa), b, e, c, p, center[0], center[2], tangens*180./np.pi, turn, np.abs(r1*r2 / ( r1 + r2 )), reflectionAngle, r1, r2
# internal methods to calculate missing data from
# input waist w1,
# output waist D2
# input waist distance to surface (D1)
# output waist distance to surface (D2)
# focal length
def _calcD2 ( w1, w2, d1, wavelength) :
d2 = np.pi**2. * w2**2. / (wavelength**2.) * ( w1**2. - w2**2.)
d2 += d1**2.* w2**2. / (w1**2.)
d2 /= np.sqrt(np.absolute(d2))
return d2
def _calcW2 ( d1, d2, w1, wavelength) :
p = - w1**2 - wavelength**2.*d1**2./(np.pi**2.*w1**2.)
q = wavelength**2.*d2**2/(np.pi**2.)
if p**2. - q < 0. :
return[None, None]
try :
w2 = [ np.sqrt(-p/2. + np.sqrt(p**2./4.-q)), np.sqrt(-p/2. - np.sqrt(p**2./4.-q))]
except :
return [None, None]
#check if one of the solutions is zero (or close to)
if w2 [0] < 0.0001 : w2[0] = w2[1]
if w2 [1] < 0.0001 : w2[1] = w2[0]
return [np.amin(w2), np.amax(w2)]
def _calcF ( w1, w2, d1, d2, wavelength) :
zr1 = np.pi * w1 * w1 / wavelength
zr2 = np.pi * w2 * w2 / wavelength
if d1 == 0. : d1 = 1e-6
if d2 == 0. : d2 = 1e-6
r1 = d1 + zr1*zr1 / d1
r2 = d2 + zr2*zr2 / d2
return r1*r2 / ( r1 + r2 )
def _calcW2D2 ( w1, d1, f, wavelength):
tt = d1/f-1.
zr = np.pi * w1 * w1 /wavelength
d2 = f * (1. + tt/(tt*tt + zr*zr/(f*f) ) )
wout = _calcW2(d1, d2, w1, wavelength)
if wout[0] == None :
return None, None
if np.absolute(_calcF(w1, wout[0], d1, d2, wavelength) - f) > 0.01 :
return wout[1], d2
return wout[0], d2
def _calcD1D2 ( w1, w2, f, wavelength, getBigW = False):
tt = -1.
if getBigW:
tt = 1.
zr = np.pi * w1 * w1 /wavelength
if ( w1*w1 / (w2*w2) - zr*zr/(f*f) < 0.) :
return None, None
d1 = f+tt*f *np.sqrt(w1*w1 / (w2*w2) - zr*zr/(f*f))
d2 = _calcD2(w1, w2, d1, wavelength)
if np.absolute(_calcF(w1, w2, d1, d2, wavelength) - f) > 0.01 :
d2 = -d2
return d1, d2
def _calcW1W2 ( d1, d2, f, wavelength):
mult = ( d1 - f )/( d2 - f ) - ( d1 / f - 1)**2.
if mult < 0 :
return None, None
w1 = sqrt(f*wavelength/np.pi * sqrt(mult))
wout = _calcW2(d1, d2, w1, wavelength)
if np.absolute(_calcF(w1, wout[0], d1, d2, wavelength) - f) > 0.01 :
return w1, wout[1]
return w1, wout[0]
def _calcW1D2 ( d1, w2, f, wavelength, getBigW = False):
zr = np.pi * w2 * w2 /wavelength
p = -f*f/(d1-f) - 2*f
q = f*f * ( f/(d1-f) + zr*zr + 1)
if (p*p/4. - q < 0) :
return None, None
d2a = - p/2.0 + sqrt( p*p / 4.0 - q )
d2b = - p/2.0 - sqrt( p*p / 4.0 - q )
win = _calcW2 ( d2, d1, w2, wavelength)
l1 = np.absolute(_calcF(w1, win[0], d1, d2a, wavelength) - f)
l2 = np.absolute(_calcF(w1, win[1], d1, d2a, wavelength) - f)
l3 = np.absolute(_calcF(w1, win[0], d1, d2b, wavelength) - f)
l4 = np.absolute(_calcF(w1, win[1], d1, d2b, wavelength) - f)
if getBigW:
if l2 < 0.01 :
return win[1], d2a
if l4 < 0.01 :
return win[1], d2b
if l1 < 0.01 :
return win[0], d2a
return win[0], d2b
else :
if l3 < 0.01 :
return win[0], d2b
if l1 < 0.01 :
return win[0], d2a
if l4 < 0.01 :
return win[1], d2b
return win[1], d2a
# method to complete data set using three of the above mentioned 5 values. Methos internally uses the above functions
# use None to indcate values to be calculated
[docs]
def complete ( wavelength, w1 = None, w2 = None, d1 = None, d2 = None, f = None, getBigW = False) :
# f not available
if ( w1 == None ) and ( w2 != None ) and ( d1 != None ) and ( d2 != None ) and ( f == None ):
w1 = _calcW2(d2, d1, w2, wavelength)
if getBigW :
return [w1[1], w2, d1, d2, _calcF(w1[1], w2, d1, d2, wavelength)]
return [w1[0], w2, d1, d2, _calcF(w1[0], w2, d1, d2, wavelength)]
if ( w2 == None ) and ( w1 != None ) and ( d1 != None ) and ( d2 != None ) and ( f == None ):
w2 = _calcW2(d1, d2, w1, wavelength)
if getBigW :
return [w1, w2[1], d1, d2, _calcF(w1, w2[1], d1, d2, wavelength)]
return [w1, w2[0], d1, d2, _calcF(w1, w2[0], d1, d2, wavelength)]
if ( d1 == None ) and ( d2 != None ) and ( w1 != None ) and ( w2 != None ) and ( f == None ):
d1 = _calcD2(w2, w1, d2, wavelength)
return [w1, w2, d1, d2, _calcF(w1, w2, d1, d2, wavelength)]
if ( d2 == None ) and ( d1 != None ) and ( w1 != None ) and ( w2 != None ) and ( f == None ):
d2 = _calcD2(w1, w2, d1, wavelength)
return [w1, w2, d1, d2, _calcF(w1, w2, d1, d2, wavelength)]
# f avalable
if ( w1 != None ) and ( d1 != None ) and ( w2 == None ) and ( d2 == None ) and ( f != None ):
w2, d2 = _calcW2D2(w1, d1, f, wavelength)
return [w1, w2, d1, d2, f]
if ( w1 == None ) and ( d1 == None ) and ( w2 != None ) and ( d2 != None ) and ( f != None ):
w1, d1 = _calcW2D2(w2, d2, f, wavelength)
return [w1, w2, d1, d2, f]
if ( w1 != None ) and ( d1 == None ) and ( w2 != None ) and ( d2 == None ) and ( f != None ):
d1, d2 = _calcD1D2(w1, w2, f, wavelength, getBigW)
return [w1, w2, d1, d2, f]
if ( w1 == None ) and ( d1 != None ) and ( w2 == None ) and ( d2 != None ) and ( f != None ):
w1, w2 = _calcW1W2(d1, d2, f, wavelength)
return [w1, w2, d1, d2, f]
if ( w1 == None ) and ( d1 != None ) and ( w2 != None ) and ( d2 == None ) and ( f != None ):
w1, d2 = _calcW1D2(d1, w2, f, wavelength, getBigW)
return [w1, w2, d1, d2, f]
if ( w1 != None ) and ( d1 == None ) and ( w2 == None ) and ( d2 != None ) and ( f != None ):
w2, d1 = _calcW1D2(d2, w1, f, wavelength, getBigW)
return [w1, w2, d1, d2, f]
return [None, None, None, None, None]