Commit 87ba9bab authored by Daniel Toyra's avatar Daniel Toyra

Updated Virgo tools

parent e58d6d94
......@@ -103,9 +103,9 @@ def round_to_n(x, n):
factor = (10 ** power)
return round(x * factor) / factor
def vprint(verbose, printstr):
def vprint(verbose, printstr, end='\n'):
if verbose:
print(printstr)
print(printstr,end=end)
def BS_optical_path(thickness, n=1.44963098985906, angle=45.0):
"""
......
......@@ -493,7 +493,7 @@ class ADV_IFO(IFO):
else:
raise pkex.BasePyKatException("\033[91m offset_type must be DARM or MICH. \033[0m")
print("-- applying user-defined DC offset to {}:".format(offset_type))
vprint(verbose, "-- applying user-defined DC offset to {}:".format(offset_type))
_kat = self.kat
m = self.mirrors
......@@ -533,11 +533,11 @@ class ADV_IFO(IFO):
out = kat.run()
print("-- adjusting {} DCoffset based on light in dark port:".format(offset_type))
vprint(verbose, "-- adjusting {} DCoffset based on light in dark port:".format(offset_type))
waste_light = round(float(out[signame]),1)
print(" waste light in AS port of {:2} W".format(waste_light))
vprint(verbose, " waste light in AS port of {:2} W".format(waste_light))
#kat_lock = _kat.deepcopy()
DCoffset = self.find_DC_offset(5*waste_light, offset_type, verbose=verbose)
......@@ -551,6 +551,7 @@ class ADV_IFO(IFO):
This function directly alters the tunings of the associated kat object.
"""
m = self.mirrors
if offset_type == 'DARM' or offset_type == 'darm':
isDARM = True
......@@ -603,7 +604,7 @@ class ADV_IFO(IFO):
#kat.NI.phi = IXphi - phi/2.0
out = kat.run()
print(" ! ", out[self.B1.get_signal_name()], phi)
# print(verbose, " ! ", out[self.B1.get_signal_name()], phi)
return np.abs(out[self.B1.get_signal_name()] - AS_power)
......@@ -881,6 +882,109 @@ class ADV_IFO(IFO):
elif k == mirrors['IY']:
self.kat.CPW_TL.f = float(v['f_CP_new'])
def find_maxtem(self, tol=5e-3, verbose=False):
'''
Finding the minimum required maxtem for the power to converge to within the relative tolerance tol.
'''
kat1 = self.kat.deepcopy()
sigs = []
sigs.append(kat1.IFO.POW_BS.add_signal())
sigs.append(kat1.IFO.POW_X.add_signal())
sigs.append(kat1.IFO.POW_Y.add_signal())
kat1.parse('noxaxis\nyaxis abs')
run = True
P_old = np.zeros([2,3], dtype=float) + 1e-9
P = np.zeros(3, dtype=float) + 1e-9
rdiff = np.ones([2,3],dtype=float)
mxtm = 0
while run:
P_old[0,:] = P_old[1,:]
P_old[1,:] = P
rdiff_old = rdiff[1,:]
kat1.maxtem = mxtm
out = kat1.run()
# print(out.stdout)
for k,s in enumerate(sigs):
P[k] = out[s]
rdiff = np.abs((P-P_old)/P_old)
if rdiff.max()<tol:
run = False
# print(kat1.maxtem, rdiff, rdiff.max())
mxtm += 1
# Stepping back to the lowest acceptable maxtem
mxtm -= 2
# One more step back if the two previous iterations were within the tolerance
if rdiff_old.max() < tol:
mxtm -= 1
self.kat.maxtem = mxtm
vprint(verbose, "Maxtem set to {}".format(mxtm))
def find_warm_detector(self, mirror_list, DCoffset, verbose=False):
"""
Computes the thermal effects for the mirrors specified in mirror_list and sets the
warm interferometer values. For an input test masse, the thermal lens is computed
and is combined with the CP-lens into a new CP-lens.
mirror_list - List of mirror names to compute the thermal effect for
DCoffset - The DARM DC-offset used in this file [degrees]
verbose - If set to True, information is printed.
"""
kat1 = self.kat.deepcopy()
# Cold IFO RoCs and focal lengthts
new = copy.copy(kat1.IFO.cold_ifo)
tol = 1e-3
run = True
a = 0
while run:
a += 1
vprint(verbose, 'Iteration {}'.format(a))
old = copy.copy(new)
# Computing new paramaters
vprint(verbose, ' Finding maxtem...',end=" ")
kat2 = kat1.deepcopy()
kat2.IFO.remove_modulators()
kat2.IFO.find_maxtem(tol=1e-3)
vprint(verbose, '{}\n Re-tuning interferometer...'.format(kat2.maxtem), end=" ")
pretune(kat2, 1e-7, verbose=False)
# adv.pretune_status(kat2)
kat1.IFO.apply_tunings(kat2.IFO.get_tunings())
kat1.maxtem = kat2.maxtem
kat1.IFO.set_DC_offset(DCoffset=DCoffset, verbose=False)
vprint(verbose, 'Done!\n Computing thermal effect...', end=" ")
new, out = compute_thermal_effect(kat1, mirror_list, nScale=True)
vprint(verbose, 'Done!')
diff = np.zeros(len(new), dtype=float)
# Updating IFO parameters
for i,(k,v) in enumerate(new.items()):
if isinstance(kat1.components[k], pykat.components.lens):
kat1.components[k].f = v
vprint(verbose, ' {}.f: {:.5e} m --> {:.5e} m'.format(k,old[k],kat1.components[k].f.value), end=",")
elif (isinstance(kat1.components[k], pykat.components.mirror) or
isinstance(kat1.components[k], pykat.components.beamSplitter)):
kat1.components[k].Rc = v
vprint(verbose, ' {}.Rc: {:.5e} m --> {:.5e} m'.format(k,old[k],kat1.components[k].Rc.value), end=",")
# Comparing the new and previous parameters
diff[i] = np.abs((new[k] - old[k])/old[k])
vprint(verbose,"")
if diff.max() < tol:
run = False
vprint(verbose, 'Converged!')
# Setting new parameters to the kat-object
for i,(k,v) in enumerate(new.items()):
if isinstance(kat1.components[k], pykat.components.lens):
self.kat.components[k].f = v
elif (isinstance(kat1.components[k], pykat.components.mirror) or
isinstance(kat1.components[k], pykat.components.beamSplitter)):
self.kat.components[k].Rc = v
def assert_adv_ifo_kat(kat):
#print(ADV_IFO)
......@@ -1229,6 +1333,16 @@ def make_kat(name="design_PR", katfile=None, verbose = False, debug=False, keepC
MPs[mirrors['IY']]['aSub'] = 3.0e-5
kat.IFO.mirror_properties = MPs
# Storing RoCs and focal lengths for the cold IFO. To be used when computing thermal effects
cold = {}
for k,v in kat.components.items():
if isinstance(v, pykat.components.mirror) or isinstance(v, pykat.components.beamSplitter):
cold[k] = v.Rc.value
elif isinstance(v, pykat.components.lens):
cold[k] = v.f.value
kat.IFO.cold_ifo = cold
kat.IFO.update()
kat.IFO.lockNames = None
......@@ -1257,9 +1371,9 @@ def pretune(_kat, pretune_precision=1.0e-4, verbose=False):
# kat and associated IFO object passed in
IFO = _kat.IFO
m = IFO.mirrors
print("-- pretuning interferometer to precision {0:2g} deg = {1:2g} m".format(pretune_precision,
pretune_precision*_kat.lambda0/360.0))
vprint(verbose,"-- pretuning interferometer to precision {0:2g} deg = {1:2g} m".format(pretune_precision,
pretune_precision*_kat.lambda0/360.0))
kat = _kat.deepcopy()
kat.removeBlock("locks", False)
......@@ -1326,7 +1440,7 @@ def pretune(_kat, pretune_precision=1.0e-4, verbose=False):
vprint(verbose, " found max/min at: {} (precision = {:2g})".format(phi, precision))
IFO.preSRCL.apply_tuning(phi)
print(" ... done")
vprint(verbose," ... done")
......@@ -1524,7 +1638,7 @@ def thermal_lensing(kat, mirror_list, nScale=False):
kat1.CPN_TL.f = float(v['f_CP_new'])
elif k == mirrors['IY']:
kat1.CPW_TL.f = float(v['f_CP_new'])
return kat1
return
def compute_thermal_effect(kat, mirror_list, nScale=False):
......@@ -1596,6 +1710,9 @@ def compute_thermal_effect(kat, mirror_list, nScale=False):
#################################
# Compute thermal effects
#################################
# Dictionary for storing new component parameters
new_params = {}
# Dictionary for storing auxiliary data
output = {}
for k,v in Ms.items():
res = {}
......@@ -1637,24 +1754,24 @@ def compute_thermal_effect(kat, mirror_list, nScale=False):
res[k] += (f,)
'''
# Combining CP and input mirror thermal lenses into 1
# Combining CP and input mirror thermal lenses into one new lens at the CP
if k == mirrors['IX']:
# Initial focal length of CP
f_old = kat1.CPN_TL.f.value
# Distance between CP and input mirror
d = kat1.sCPN_NI.L.value
new_params['CPN_TL'], errors = combine(kat1.IFO.cold_ifo['CPN_TL'], res['f_thermal'], d=d, q=q)
elif k == mirrors['IY']:
# Initial focal length of CP
f_old = kat1.CPW_TL.f.value
# Distance between CP and input mirror
d = kat1.sCPW_WI.L.value
# Combining CP and input mirror lenses into one new lens at the CP
f_CP_new, errors = combine(f_old, res['f_thermal'], d=d, q=q)
res['f_CP_new'] = f_CP_new
res['rel_errs'] = errors
new_params['CPW_TL'], errors = combine(kat1.IFO.cold_ifo['CPW_TL'], res['f_thermal'], d=d, q=q)
res['compound_lens_errs'] = errors
output[k] = res
return output
return new_params, output
......@@ -12,9 +12,9 @@ dtoyra@star.sr.bham.ac.uk
# Notes to self
# - Several references missing for claimed measured values.
# - Extract all thermal lenses and make them into one thermal lens at
# each compensation plate.
# - Tools to make several lenses into one
# - Function computing thermal lens from arm cavity power, absorption
# coefficients and Roc of CP surfaces
# - Python scirpt computing Finesse and rt-loss of a pretuned kat-file.
......@@ -22,11 +22,34 @@ dtoyra@star.sr.bham.ac.uk
# or carrier gain (not sure that works in general)
# - Prepare file that computes arm cavity finesse as a function of etalon
# tunings.
# - Several references missing for claimed measured values.
# - After having thermal lens from Hello-Vinet, Plot PR-gain over thermal lens offset, and over PR-ROC.
# - s
# Setup file thermal effects file with maxmimised power at 14 W
# -------------------------------------------------------------
# 1. Set input power to 14 W.
# 2. Compute thermal effects given the powers for 14 W input power,
# for the well tuned and matched file.
# 3. Set new RoCs and IM thermal lenses of the kat-object according to the thermal
# effects obtained.
# 4. "Mode match"/maximise powers by tuning AR_RoC of IMs and CP lenses.
# 5. Store AR_RoCs and CP focal lengths in a permanent kat-file.
# 6. Store warm parameters in adv-tools file.
# Thermal effect computation (kat-file assumed to be prepared properly):
# ----------------------------------------------------------------------
# 1. Chose an input power (0 to 14 W should be tested. Should bring the IFO to the
# state obtained above. For other values, start from 14 W and move towards
# new power)
# 2. Compute new RoCs and thermal lenses, and set these to kat-object.
# 3. Tune file, check maxtem needed, and obtain new powers.
# 4. Repeat 2 and 3 until tolerance reached.
%%% FTblock Laser
###########################################################################
# The input power used during last Science run
l i1 14 0 0 nin
s s0 1m nin nEOM1a
......@@ -187,7 +210,7 @@ s LW 2999.8 nWI2 nWE1
# AR-reflectivity is set as a loss.
m1 WE $T_WE $L_WE $WE_phi nWE1 nWEsub1
s sWEsub .2 $nsilica nWEsub1 nWEsub2
m2 WEAR 0 $L_WEAR $WEAR_phi nWEsub2 nWE2
m2 WEAR 0 $L_WEAR 0 nWEsub2 nWE2
###########################################################################
%%% FTend Warm
......@@ -202,7 +225,7 @@ m2 WEAR 0 $L_WEAR $WEAR_phi nWEsub2 nWE2
s lsr 5.943 nBSe nSR1
m SR 0 1 0 nSR1 nSRsub1
s SRsub 0.1 $nsilica nSRsub1 nSRsub2
s sSRsub 0.1 $nsilica nSRsub1 nSRsub2
m SRAR 0 1 0 nSRsub2 nSR2
###########################################################################
......@@ -290,12 +313,13 @@ s sOMC2_OC_CS2 0.0630 $nsilica nOMC2_OCb nOMC2_CS2a
bs OMC2_CS2 0.9999665 33.5u 0 8.876 nOMC2_CS2a nOMC2_CS2b nOMC2_CS2c nOMC2_CS2d
s sOMC2_CS2_IC 0.0600 $nsilica nOMC2_CS2b nOMC2_ICd
# Space to output B1
s sOMC2_B1 0 nOMC2_OCc nB1
###########################################################################
%%% FTend OMC
# Space to output B1
s sOut 0 nOMC2_OCc nB1
%%% FTblock Gaussian
###########################################################################
......@@ -338,10 +362,7 @@ attr OMC2_CS2 Rc 1.499
# -----------
# attr WIAR Rc -1420 # Designed WI AR RoC [TDR, table 2.6]
# attr NIAR Rc -1420 # Designed NI AR RoC [TDR, table 2.6]
attr WIAR Rc -1423.95 # Approximately the same as measured HR surface,
# and slightly adjusted so that a mode matched file
# model can be obtained with equal CP thermal
# lensing.
attr WIAR Rc -1424.5 # Approximately the same as measured HR surface
attr NIAR Rc -1424.6 # Approximately the same as measured HR surface
attr PRAR Rc -3.62 # Measured PR AR RoC [VIR-0029A-15]
attr SRAR Rc 3.59 # Design [TDR, table 2.8]
......@@ -349,8 +370,12 @@ attr SRAR Rc 3.59 # Design [TDR, table 2.8]
# Lenses
# -----------
# Compensation plate thermal lenses
const f_CPN_TL 50890 # North
const f_CPW_TL 50890 # West
const f_CPN_TL 60509.52414597 # North
const f_CPW_TL 64643.60054991 # West
#const f_CPN_TL 54000 # North
#const f_CPW_TL 930000 # West
###########################################################################
%%% FTend ROCs
......
......@@ -18,6 +18,8 @@ import scipy.optimize as so
import pylab
from scipy.optimize import minimize
import pykat.exceptions as pkex
from pykat.tools.lensmaker import lensmaker
def hellovinet(P_coat, P_sub_in, P_sub_out, HR_RoC = None, AR_RoC = None, HR_zOff=None, AR_zOff=None, **kwargs):
......@@ -59,6 +61,9 @@ def hellovinet(P_coat, P_sub_in, P_sub_out, HR_RoC = None, AR_RoC = None, HR_zOf
P_coat - Power on HR coating from vacuum side. See figure above.
P_sub_in - Power propagaint from AR-surface to HR-surface. See figure above.
P_sub_out - Power propagating from HR-surface to AR-surface. See figure above.
mirror_properties - Dictionary where all the below properties can be specified.
HR_RoC - Initial RoC of the HR surface [m]. If specified, the new RoC will be
fitted to the deformed surface. For sign, see above.
AR_RoC - Same as for HR_RoC, but for the AR-surface.
......@@ -67,9 +72,6 @@ def hellovinet(P_coat, P_sub_in, P_sub_out, HR_RoC = None, AR_RoC = None, HR_zOf
which affects the RoC fitting. Set this to None (no fitting) or 0 (the
offset should always be small).
AR_zOff - Same as for HR_zOff above, but for the AR-surface.
mirror_properties - Dictionary where all the below properties can be specified directly.
thickness - Mirror thickness [m] (default 0.2 m)
aCoat - Coating power absorption (default 2.5 ppm)
aSub - Substrate power absorptio [1/m] default (30 ppm/m)
......@@ -84,6 +86,7 @@ def hellovinet(P_coat, P_sub_in, P_sub_out, HR_RoC = None, AR_RoC = None, HR_zOf
dndT - Mirror index of refraction change with temperature (default 8.7 ppm)
N - Number of data points along the mirror radius a (default 176).
nScale - If set to true, the optical path length data is scaled to physical distance.
fitCurv - If set to true, the curvature is fitted instead of the RoC.
Returns:
--------
......@@ -106,28 +109,28 @@ def hellovinet(P_coat, P_sub_in, P_sub_out, HR_RoC = None, AR_RoC = None, HR_zOf
#############################################
sigmab = 5.670367e-8 # Stephan-Boltzmann constant [W m**(-2) K**(-4)]
h=0.2 # test mass thickness [m]
aCoat=2.5e-6 # coating absorption
aSub=3.0e-5 # substrate absorption [1/m] [TDR, table 2.6]
n=1.452 # SiO2 refraction index
a=0.175 # test mass radius [m]
w=0.049 # Spot size at Virgo input mirrors
K=1.380 # SiO2 thermal conductivity
T0=295.0 # room temperature
emiss=0.89 # SiO2 emissivity
alpha=0.54e-6 # SiO2 thermal expansion coefficient
sigma=0.164 # SiO2 Poisson's ratio
dndT=8.7e-6 # Mirror index of refraction change with temperature (default 8.7 ppm)
N = 176 # Number of data points along the mirror radius
scale = 1.0 # Scales path-length data before fitting RoCs
nScale = False # If true, scale is set to 1/n
h=0.2 # test mass thickness [m]
aCoat=2.5e-6 # coating absorption
aSub=3.0e-5 # substrate absorption [1/m] [TDR, table 2.6]
n=1.452 # SiO2 refraction index
a=0.175 # test mass radius [m]
w=0.049 # Spot size at Virgo input mirrors
K=1.380 # SiO2 thermal conductivity
T0=295.0 # room temperature
emiss=0.89 # SiO2 emissivity
alpha=0.54e-6 # SiO2 thermal expansion coefficient
sigma=0.164 # SiO2 Poisson's ratio
dndT=8.7e-6 # Mirror index of refraction change with temperature (default 8.7 ppm)
N = 176 # Number of data points along the mirror radius
scale = 1.0 # Scales path-length data before fitting RoCs
nScale = False # If true, scale is set to 1/n
fitCurv = False # If ture, the curvature is fitted instead of the RoC.
#Pin=125; # IFO input power
#RecG=37.5; # recycling gain
#Finesse=443; # F-P cavity finesse not used
# Updating values specified as a dictionary.
if 'mirror_properties' in kwargs:
for k,v in kwargs['mirror_properties'].items():
......@@ -159,7 +162,17 @@ def hellovinet(P_coat, P_sub_in, P_sub_out, HR_RoC = None, AR_RoC = None, HR_zOf
N = v
elif k == 'nScale':
nScale = v
elif k == 'HR_RoC':
HR_RoC = v
elif k == 'AR_RoC':
AR_RoC = v
elif k == 'HR_zOff':
HR_zOff = v
elif k == 'AR_zOff':
AR_zOff = v
elif k == 'fitCurv':
fitCurv = v
# Updating values specified as arguments. These have precedence
# over the parameters specified in the dictionary.
for k, v in kwargs.items():
......@@ -191,9 +204,9 @@ def hellovinet(P_coat, P_sub_in, P_sub_out, HR_RoC = None, AR_RoC = None, HR_zOf
N = v
elif k == 'nScale':
nScale = v
elif k == 'fitCurv':
fitCurv = v
# Scales the optical path length data to physical distances
if nScale:
scale = 1.0/n
......@@ -266,7 +279,7 @@ def hellovinet(P_coat, P_sub_in, P_sub_out, HR_RoC = None, AR_RoC = None, HR_zOf
######################################################################################
# To Valeria:
# ----------
# I take responsibility for everything that below this line, so you don't
# I take responsibility for everything below this line, so you don't
# need to check anything below here. /Daniel
######################################################################################
......@@ -280,42 +293,57 @@ def hellovinet(P_coat, P_sub_in, P_sub_out, HR_RoC = None, AR_RoC = None, HR_zOf
# should affect 1. The other three contributions gets added as a
# thermal lens.
# Setting HR-deformation
OPL_HR = OPLtec*scale
# Setting AR-thermal
OPL_AR = (OPLs + OPLtes + OPLc)*scale
fitData = []
# Computing new RoC for HR-surface, if initial HR_RoC was given.
if not HR_RoC is None:
# Creating initial mirror surface
Z_HR0 = createSurface(r, HR_RoC, HR_zOff)
# Adding the distortion
Z_HR1 = Z_HR0 + OPL_HR
# Fitting
out = fit_circle(r, Z_HR1, Rc0 = HR_RoC, zOff0 = HR_zOff, w = w)
fitData.append(out)
#if isinstance(out, list):
# HR_RoC1 = out[0]
# HR_zOff1 = out[1]
#else:
# HR_RoC1 = out
# Setting HR-deformation (not used)
# OPL_HR = OPLtec*scale
# Setting AR-thermal deformation
OPL_AR = (OPLc + OPLtec + OPLs + OPLtes)*scale
# Thickness of the equivalent lens
d = OPL_AR[np.abs(r).argmin()]
# Computing new RoC for HR-surface, if initial HR_RoC was given. (not used)
#if not HR_RoC is None:
# # Creating initial mirror surface
# Z_HR0 = createSurface(r, HR_RoC, HR_zOff)
# # Adding the distortion
# Z_HR1 = Z_HR0 + OPL_HR
# # Fitting
# out = fit_circle(r, Z_HR1, Rc0 = HR_RoC, zOff0 = HR_zOff, w = w)
# rocData.append(out)
# #if isinstance(out, list):
# # HR_RoC1 = out[0]
# # HR_zOff1 = out[1]
# #else:
# # HR_RoC1 = out
# Computing new effective RoC for AR-surface, if initial AR_RoC was given.
if not AR_RoC is None:
######################
# AR-surface
######################
# Creating initial mirror surface
Z_AR0 = createSurface(r, AR_RoC, AR_zOff)
# Adding the distortion
Z_AR1 = Z_AR0 - OPL_AR
# Fitting
out = fit_circle(r, Z_AR1, Rc0 = AR_RoC, zOff0 = AR_zOff, w = w)
fitData.append(out)
#if not AR_RoC is None:
# ######################
# # AR-surface
# ######################
# # Creating initial mirror surface
# Z_AR0 = createSurface(r, AR_RoC, AR_zOff)
# # Adding the distortion
# Z_AR1 = Z_AR0 - OPL_AR
# # Fitting
# rc = fit_circle(r, Z_AR1, Rc0 = AR_RoC, zOff0 = AR_zOff, w = w)
# rocData.append(rc)
# f = lensmaker(rc, AR_RoC, d, n)
# lensData.append(f)
# Computing curvature and focal length of the effective curved-flat lens.
if fitCurv:
c = fit_curvature(r, OPL_AR, Rc0 = AR_RoC, w = w)
rc = 1.0/c
else:
out = fit_circle(r, OPL_AR, Rc0 = -2000, zOff0 = AR_zOff, w = w)
rc = out[0]
return r, fitData, oplData
f = lensmaker(np.abs(rc), np.inf, d, n)
return f, [r, oplData, rc, d]
def createSurface(r,Rc,zOffset=None):
......@@ -423,6 +451,71 @@ def fit_circle(r, z, Rc0=None, zOff0=None, w=None, maxfev=2000):
if not out['success']:
msg = ' Warning: ' + out['message'].split('.')[0] + ' (nfev={:d}).'.format(out['nfev'])
return out['x']
return out.x
def fit_curvature(r, z, Rc0=None, w=None):
'''
Fits a circle to the data (r,z) minimising the difference in curvature (1/RoC) between
the cicle and the data (r, z).
r - Array with radial distances away from from the center of the mirror [m].
z - Array with mirror heights (along optical axis) at the positions r [m].
Rc0 - Initial guess of curvature [1/m].
w - Gaussian weights parameter [m]. Should be equal to beam spot radius.
'''
N = len(z)
dx = r[1]-r[0]
dz = np.zeros(N)
ddz = np.zeros(N)
# First derivatie with central difference of 4th order accuracy
dz[2:-2] = (z[:-4]/4.0 - 2.0*z[1:-3] + 2.0*z[3:-1] - z[4:]/4.0)/(3.0*dx)
# First derivative with central difference of 2nd order accuracy
dz[1] = (z[2] - z[0])/(2.0*dx)
dz[-2] = (z[-1] - z[-3])/(2.0*dx)
# First derivative with forward difference of 2nd order accuracy
dz[0] = (-3.0*z[0] + 4.0*z[1] - z[2])/(2.0*dx)
# First derivative with backward difference of 2nd order accuracy
dz[-1] = (3.0*z[-1] - 4.0*z[-2] + z[-3])/(2.0*dx)
# Second derivative with central difference of 4h order accuracy
ddz[2:-2] = (-z[:-4]/4.0 + 4.0*z[1:-3] - (15.0/2.0)*z[2:-2] + 4.0*z[3:-1] - z[4:]/4.0)/(3.0*dx**2)
# Second derivative with central difference of 2h order accuracy
ddz[1] = (z[0] - 2.0*z[1] + z[2])/dx**2
ddz[-2] = (z[-3] - 2.0*z[-2] + z[-1])/dx**2
# Second derivative with forward difference of 2nd order accuracy
ddz[0] = (2.0*z[0] - 5.0*z[1] + 4.0*z[2] - z[3] )/dx**2
# Second derivative with backward difference of 2nd order accuracy
ddz[-1] = (2.0*z[-1] - 5.0*z[-2] + 4.0*z[-3] - z[-4] )/dx**2
# Computing curvature
curv = ddz/(1.0+dz**2)**(1.5)
# Setting weights
if w is None:
weight = 1.0
else:
weight = np.exp(-2*r**2/w**2)
def func(c):
return (weight*(curv-c)**2).sum()/np.abs(weight*curv).sum()
# Initial guess of best fitting curvature
if not Rc0 is None:
p = (Rc0,)
else:
p = (0,)
opts = {'xtol': 1.0e-9, 'ftol': 1.0e-9, 'maxiter': 10000, 'maxfev': 2000, 'disp': False}
out = minimize(func, p, method='Nelder-Mead', options=opts)
if not out.success:
pkex.printWarning(out.message)
return out.x[0]
# return dz, ddz, curv, out