|
new file 100644
|
|
|
#! /cwi/bin/python
|
|
|
|
|
|
from csplot import *
|
|
|
import sys
|
|
|
|
|
|
# Update and use the corresponding SciPy functions
|
|
|
#from Scientific.Functions import Interpolation
|
|
|
import math
|
|
|
|
|
|
import numpy as num
|
|
|
import scipy as sci
|
|
|
from scipy.optimize.optimize import fmin_powell as fmin
|
|
|
from scipy.optimize.minpack import leastsq
|
|
|
from optparse import OptionParser, Option, OptionValueError
|
|
|
|
|
|
INFTY = 1e10
|
|
|
|
|
|
try:
|
|
|
BASE_PATH = os.environ['CSTREAM_DEFAULT_OUT_PATH']
|
|
|
except KeyError:
|
|
|
BASE_PATH = '/export/scratch1/luque/cstream/'
|
|
|
|
|
|
def max_in_axis(lz, m, sign, zmin=None, zmax=None):
|
|
|
""" Calculates and returns the maximum of a var in the axis and
|
|
|
its position. """
|
|
|
|
|
|
max_i = -1
|
|
|
|
|
|
# No infinity symbol in some old versions of Python
|
|
|
max_v = 0
|
|
|
min_v = 1e100
|
|
|
corr = 0
|
|
|
|
|
|
dz = lz[1] - lz[0]
|
|
|
|
|
|
if zmin is None:
|
|
|
zmin = 0
|
|
|
|
|
|
if zmax is None:
|
|
|
zmax = INFTY
|
|
|
|
|
|
for i in range(len(lz)):
|
|
|
if lz[i] < zmin or lz[i] > zmax:
|
|
|
continue
|
|
|
|
|
|
v = m[i]
|
|
|
if abs(v) > abs(max_v) and sign * v >= 0:
|
|
|
max_i = i
|
|
|
max_v = v
|
|
|
|
|
|
if abs(v) < abs(min_v) and sign * v >= 0:
|
|
|
min_i = i
|
|
|
min_v = v
|
|
|
|
|
|
if max_i >= 1:
|
|
|
corr, max_v = quadratic_max(m[max_i - 1: max_i + 2])
|
|
|
|
|
|
return lz[max_i] + corr * dz, max_v, min_v
|
|
|
|
|
|
|
|
|
def max_off_axis(lr, lz, m, sign):
|
|
|
Z, R = num.meshgrid(lz, lr)
|
|
|
indx = (sign * m).argmax(axis=None)
|
|
|
return R.flat[indx], Z.flat[indx]
|
|
|
|
|
|
|
|
|
def find_front(var, grid, path=BASE_PATH, z0=None, ret_z0=False,
|
|
|
fmin=0, fmax=INFTY, rmin=None, rmax=INFTY,
|
|
|
pre_margin=400, post_margin=20, zmin=None, zmax=None,
|
|
|
shift_origin=True, sign=None):
|
|
|
|
|
|
|
|
|
if sign is None:
|
|
|
dict_signs = {'charge': -1, 'ez': 0, 'eabs': 0}
|
|
|
sign = dict_signs.get(var, 0)
|
|
|
|
|
|
data = []
|
|
|
lr, lz, m = load_var(var, grid, path=path)
|
|
|
|
|
|
if rmin is None:
|
|
|
rmin = 0.5 * (lr[1] + lr[2])
|
|
|
|
|
|
if z0 is None:
|
|
|
z0, max_var, _ = max_in_axis(lz, m[0, :], sign,
|
|
|
zmin=zmin, zmax=zmax)
|
|
|
|
|
|
fmin, fmax = abs(max_var) * fmin, abs(max_var) * fmax
|
|
|
print fmin, fmax
|
|
|
|
|
|
|
|
|
if shift_origin:
|
|
|
z0_shift = z0
|
|
|
else:
|
|
|
z0_shift = 0
|
|
|
|
|
|
for j in xrange(len(lz)):
|
|
|
# This is a small WTF, really
|
|
|
if lz[j] < z0 - pre_margin or lz[j] > z0 + post_margin:
|
|
|
continue
|
|
|
|
|
|
r, f = find_max_row(lr, sign * m[:, j])
|
|
|
|
|
|
if abs(f) >= fmin and abs(f) < fmax and r > rmin and r < rmax:
|
|
|
data.append((lz[j] - z0_shift, r, f))
|
|
|
#else:
|
|
|
# print ("f out of (%f, %f) or r out of (%f, %f)" %
|
|
|
# (fmin, fmax, rmin, rmax))
|
|
|
|
|
|
if ret_z0:
|
|
|
return data, z0
|
|
|
else:
|
|
|
return data
|
|
|
|
|
|
class ThresholdNotFound(ValueError):
|
|
|
pass
|
|
|
|
|
|
def find_envelope(grid, path=BASE_PATH, eps=0.0005, tail=350.0):
|
|
|
var = 'charge'
|
|
|
lr, lz, m = [num.array(a) for a in load_var(var, grid, path=path)]
|
|
|
|
|
|
m = -m
|
|
|
front = num.zeros((len(lz), 2), 'd')
|
|
|
found = num.array([False for i in xrange(len(lz))])
|
|
|
|
|
|
dz = lz[1] - lz[0]
|
|
|
ntail = int(tail / dz)
|
|
|
|
|
|
lr_inv = lr[-1::-1]
|
|
|
|
|
|
for i in xrange(len(lz)):
|
|
|
try:
|
|
|
a, rmax = find_threshold(eps, m[-1::-1, i], lr_inv)
|
|
|
front[i, 0:2] = lz[i], rmax
|
|
|
found[i] = True
|
|
|
except ThresholdNotFound:
|
|
|
pass
|
|
|
|
|
|
front = front.compress(found, axis=0)
|
|
|
zlen, _ = front.shape
|
|
|
if zlen > ntail:
|
|
|
zlen = ntail
|
|
|
|
|
|
return front[-zlen:, :]
|
|
|
|
|
|
def find_threshold(threshold, a, r):
|
|
|
for ia, ir in zip(a, r):
|
|
|
if ia >= threshold:
|
|
|
return ia, ir
|
|
|
|
|
|
raise ThresholdNotFound
|
|
|
|
|
|
def simple_find_front(var, grid, path=BASE_PATH, tail=350.0, eps=1e-4):
|
|
|
var = 'charge'
|
|
|
lr, lz, m = [num.array(a) for a in load_var(var, grid, path=path)]
|
|
|
m = num.abs(m)
|
|
|
|
|
|
maxes = m.argmax(0).squeeze()
|
|
|
front = num.array([(lz[i], lr[maxes[i]])
|
|
|
for i in xrange(len(lz))
|
|
|
if maxes[i] > 0 and m[maxes[1], i] > eps])
|
|
|
|
|
|
zlen, _ = front.shape
|
|
|
if zlen > ntail:
|
|
|
zlen = ntail
|
|
|
|
|
|
return front[-zlen:, :]
|
|
|
|
|
|
|
|
|
def fit_radius(var, grid, sign, path=BASE_PATH, zmax=None):
|
|
|
front, z0 = find_front(var, grid, path=path, pre_margin=400.0, sign=sign,
|
|
|
post_margin=400.0, ret_z0=True,
|
|
|
fmax=1.0, fmin=0.5, zmax=zmax)
|
|
|
|
|
|
if not front:
|
|
|
return num.NAN
|
|
|
|
|
|
front = num.array(front)
|
|
|
|
|
|
def errors(R):
|
|
|
return num.sqrt(front[:, 1]**2 + (front[:, 0] - R)**2) - abs(R)
|
|
|
|
|
|
try:
|
|
|
x, _ = leastsq(errors, 20)
|
|
|
except TypeError:
|
|
|
raise
|
|
|
#ipshell()
|
|
|
|
|
|
return x
|
|
|
|
|
|
|
|
|
def find_max_row(lr, row):
|
|
|
max_f = 0
|
|
|
min_f = 1e100
|
|
|
max_i = 0
|
|
|
|
|
|
dr = lr[1] - lr[0]
|
|
|
|
|
|
for i in xrange(len(lr)):
|
|
|
f = row[i]
|
|
|
if abs(f) > abs(max_f):
|
|
|
max_i = i
|
|
|
max_f = f
|
|
|
|
|
|
if max_i >= 1 and max_i < len(lr) - 1:
|
|
|
corr, max_f = quadratic_max(row[max_i - 1: max_i + 2])
|
|
|
else:
|
|
|
corr = 0
|
|
|
|
|
|
return lr[max_i] + corr * dr, max_f
|
|
|
|
|
|
|
|
|
|
|
|
def inward_norm(nx, ny, y):
|
|
|
""" gets a normal vector from the array we created making sure that
|
|
|
it points always inward"""
|
|
|
|
|
|
if num.sign(y) == num.sign(ny):
|
|
|
nx, ny = -nx, -ny
|
|
|
|
|
|
return nx, ny
|
|
|
|
|
|
def normal_to_curve(m):
|
|
|
""" given an array as [(x0, y1), ...(xN, yN)] returns
|
|
|
[(nx0, ny0), ... (nxN, nyN)] of components of a vector normal to the
|
|
|
curve """
|
|
|
print m.shape
|
|
|
data = zeros(m.shape, typecode='fd')
|
|
|
n = m.shape[0]
|
|
|
|
|
|
for i in xrange(1, n - 1):
|
|
|
x, y = m[i, 0], m[i, 1]
|
|
|
|
|
|
deltas = [(m[j, 0] - m[i, 0], m[j, 1] - m[i, 1])
|
|
|
for j in (i - 1, i + 1)]
|
|
|
nx, ny = normalize(normal(*deltas))
|
|
|
nx, ny = inward_norm(nx, ny, y)
|
|
|
|
|
|
data[i, 0] = nx
|
|
|
data[i, 1] = ny
|
|
|
|
|
|
nx, ny = normalize((-(m[1, 0] - m[0, 0]),
|
|
|
m[1, 1] - m[0, 1]))
|
|
|
|
|
|
nx, ny = inward_norm(nx, ny, m[0, 1])
|
|
|
|
|
|
data[0, 0] = nx
|
|
|
data[0, 1] = ny
|
|
|
|
|
|
nx, ny = normalize((m[n - 2, 1] - m[n - 1, 1],
|
|
|
-(m[n - 2, 0] - m[n - 1, 0])))
|
|
|
|
|
|
nx, ny = inward_norm(nx, ny, m[n - 1, 1])
|
|
|
data[n - 1, 0] = nx
|
|
|
data[n - 1, 1] = ny
|
|
|
|
|
|
print data
|
|
|
return data
|
|
|
|
|
|
def normal((x0, y0), (x1, y1)):
|
|
|
print x0, y0, x1, y1
|
|
|
|
|
|
if y0 == 0 and y1 == 0:
|
|
|
return (0.0, 1.0)
|
|
|
|
|
|
return ((-(x1**2 * y0) + (x0**2 + y0 * (y0 - y1)) * y1)
|
|
|
/ (x1 * y0 - x0 * y1),
|
|
|
(x0**2 * x1 + x1 * y0**2 - x0 * (x1**2 + y1**2))
|
|
|
/ (-(x1 * y0) + x0 * y1))
|
|
|
|
|
|
|
|
|
def normalize((x, y)):
|
|
|
a = sqrt(x**2 + y**2)
|
|
|
return (x / a, y / a)
|
|
|
|
|
|
|
|
|
def quadratic_max(y):
|
|
|
""" Finds a parabola defined by the points (dx i, y[i]) {i = 1, 2, 3}
|
|
|
and returns its maximum. """
|
|
|
|
|
|
try:
|
|
|
f1 = (y[2] - 2 * y[1] + y[0])
|
|
|
f2 = (y[2] - y[0])
|
|
|
except IndexError:
|
|
|
return num.NAN, num.NAN
|
|
|
|
|
|
if f1 == 0:
|
|
|
# The points are aligned: no parabola fit
|
|
|
return 0, max(*y)
|
|
|
|
|
|
return (- f2 / (2 * f1)), y[1] - f2**2 / f1 / 8.0
|
|
|
|
|
|
|
|
|
def var_interpolate(var, grid, r_staggered=False, z_staggered=False, **kwargs):
|
|
|
lr, lz, var = load_var(var, grid, **kwargs)
|
|
|
|
|
|
lr = num.array(lr)
|
|
|
lz = num.array(lz)
|
|
|
|
|
|
if r_staggered:
|
|
|
lr = lr + 0.5 * (lr[1] - lr[0])
|
|
|
|
|
|
if z_staggered:
|
|
|
lz = lz + 0.5 * (lz[1] - lz[0])
|
|
|
|
|
|
var = num.array(var)
|
|
|
|
|
|
lr, var = mirror_var(lr, var)
|
|
|
|
|
|
f = Interpolation.InterpolatingFunction([lr, lz], var, default=1e20)
|
|
|
return f
|
|
|
|
|
|
def find_front_iter(var, grid, n, **kwargs):
|
|
|
""" After finding a front, improves it by n iterations. """
|
|
|
print kwargs
|
|
|
|
|
|
m = find_front(var, grid, **kwargs)
|
|
|
m = num.array(m)
|
|
|
|
|
|
interp = var_interpolate(var, grid, path=kwargs.get('path', None))
|
|
|
|
|
|
for i in xrange(n):
|
|
|
norm = normal_to_curve(m)
|
|
|
m = iter_front(m, norm, interp)
|
|
|
|
|
|
return m
|
|
|
|
|
|
|
|
|
def iter_front(m, norm, interp, max_step=80, ds=0.05):
|
|
|
mnext = zeros(m.shape, typecode='fd')
|
|
|
|
|
|
for p in xrange(m.shape[0]):
|
|
|
z0, r0 = m[p, 0], m[p, 1]
|
|
|
nz, nr = norm[p, 0], norm[p, 1]
|
|
|
|
|
|
zs = [z0 + ds * i * nz for i in xrange(-max_step, max_step)]
|
|
|
rs = [r0 + ds * i * nr for i in xrange(-max_step, max_step)]
|
|
|
|
|
|
interp_zr = [interp(r, z) for r, z in zip(rs, zs)]
|
|
|
|
|
|
max_i, max_e = -1, 0
|
|
|
|
|
|
for i, e in enumerate(interp_zr):
|
|
|
print i, rs[i], zs[i], e
|
|
|
if e > max_e or max_i < 0:
|
|
|
max_i = i
|
|
|
max_e = e
|
|
|
|
|
|
mnext[p, 0] = zs[max_i]
|
|
|
mnext[p, 1] = rs[max_i]
|
|
|
|
|
|
return mnext
|
|
|
|
|
|
def fit_saffman(id, grid, var='eabs', width=64.0, fix_vertex=False, rmax=300):
|
|
|
path = os.path.join(BASE_PATH, id)
|
|
|
|
|
|
front = find_front(var, grid,
|
|
|
path=path, fmin=0.0001, rmin=1.1, rmax=rmax,
|
|
|
zmin=0, zmax=4000, shift_origin=False,
|
|
|
pre_margin=100, post_margin=5)
|
|
|
front = num.array(front)
|
|
|
|
|
|
#front = find_envelope(grid, path=path)
|
|
|
#print front
|
|
|
|
|
|
vert_x, vert_y, _ = front[num.abs(front).argmax(0)[0]]
|
|
|
#vert_x, vert_y = front[num.abs(front).argmax(0)[0]]
|
|
|
|
|
|
print "Vertex found at (%f, %f)" % (vert_x, vert_y)
|
|
|
|
|
|
L = 2 * width
|
|
|
|
|
|
def chi2((lam, C)):
|
|
|
#A = 2 * width / math.pi
|
|
|
#y = A * num.arccos(num.exp(((front[:, 0] - C) / A)).clip(-10000.0, 1.0))
|
|
|
x = front[:, 0]
|
|
|
y = (lam * L / (2 * math.pi)
|
|
|
* num.arccos((2 * num.exp(2 * math.pi * x
|
|
|
/ (L * (1 - lam)) - C) - 1.0)
|
|
|
.clip(-1.0, 1.0)))
|
|
|
|
|
|
return ((y - front[:, 1]) ** 2).sum()
|
|
|
|
|
|
if not fix_vertex:
|
|
|
lam, C = fmin(chi2, (0.5, vert_x), maxfun=100)
|
|
|
else:
|
|
|
def chi2_partial(lam):
|
|
|
return chi2((lam, vert_x))
|
|
|
|
|
|
lam, C = fmin(chi2_partial, (0.5,), maxfun=100), vert_x
|
|
|
|
|
|
|
|
|
print "\tlambda = %f\n\tx0 = %f" % (lam, C)
|
|
|
return lam, C
|
|
|
|
|
|
|
|
|
def find_emax(rid, grid, zinterval=None, path=BASE_PATH):
|
|
|
var = rid + ':eabs'
|
|
|
lr, lz, m = [num.array(a) for a in load_var(var, grid, path=path)]
|
|
|
|
|
|
z0 = lz[0]
|
|
|
dz = lz[1] - z0
|
|
|
|
|
|
if zinterval is not None:
|
|
|
zmin, zmax = zinterval
|
|
|
izmin, izmax = [int((z - z0) / dz) for z in zmin, zmax]
|
|
|
|
|
|
izmin = max(0, izmin)
|
|
|
izmax = min(m.shape[1] - 1, izmax)
|
|
|
|
|
|
m = m[:, izmin:izmax]
|
|
|
|
|
|
return max(m.flat)
|
|
|
|
|
|
def find_diameters(var, grid, neg_z, pos_z, path=BASE_PATH):
|
|
|
lr, lz, m = load_var(var, grid, path=path)
|
|
|
|
|
|
m = m * m
|
|
|
diameters = num.dot(lr, m) / num.dot(ones(lr.shape), m)
|
|
|
table = num.concatenate((lz.reshape(len(lz), 1),
|
|
|
diameters.reshape(len(lz), 1)), axis=1)
|
|
|
|
|
|
table = table.compress(num.logical_and(num.greater(lz, pos_z),
|
|
|
num.less(lz, neg_z)), axis=0)
|
|
|
#ipshell()
|
|
|
dfile = os.path.join(path, 'diameters.%s.tsv' % grid)
|
|
|
try:
|
|
|
neg_diam = num.max(table[-40:, 1])
|
|
|
except ValueError:
|
|
|
neg_diam = num.NAN
|
|
|
|
|
|
try:
|
|
|
pos_diam = num.max(table[:40, 1])
|
|
|
except ValueError:
|
|
|
pos_diam = num.NAN
|
|
|
|
|
|
num.savetxt(os.path.join(path, 'diameters.%s.tsv' % grid), table)
|
|
|
|
|
|
return neg_diam, pos_diam
|
|
|
|
|
|
def posneg_fronts(grid, path=BASE_PATH, var='charge', zmax=None):
|
|
|
lz, m = load_axis(var, grid, path=path)
|
|
|
neg_z, neg_max, neg_min = max_in_axis(lz, m, -1, zmax=zmax)
|
|
|
pos_z, pos_max, pos_min = max_in_axis(lz, m, 1, zmax=zmax)
|
|
|
return [neg_z, pos_z, neg_max, pos_max]
|
|
|
|
|
|
|
|
|
def integrate(grid, path=BASE_PATH, var='charge', sign=1.0, threshold=None,
|
|
|
zrange=None):
|
|
|
lr, lz, m = load_var(var, grid, path=path)
|
|
|
if zrange is not None:
|
|
|
z0, z1 = zrange
|
|
|
if not num.isnan(z0) and not num.isnan(z1):
|
|
|
lzmask = num.less(lz, z1) * num.greater(lz, z0)
|
|
|
m = m.compress(lzmask, axis=1)
|
|
|
lz = lz.compress(lzmask)
|
|
|
|
|
|
if threshold is not None:
|
|
|
if sign > 0:
|
|
|
thres = num.greater
|
|
|
mmax = max(num.ravel(m))
|
|
|
else:
|
|
|
thres = num.less
|
|
|
mmax = min(num.ravel(m))
|
|
|
|
|
|
m = m * thres(m, threshold * mmax)
|
|
|
|
|
|
dr = lr[1] - lr[0]
|
|
|
dz = lz[1] - lz[0]
|
|
|
|
|
|
dV = lr.reshape((lr.size, 1)) * dr * dz
|
|
|
m = m * dV
|
|
|
return sum(m.flat)
|
|
|
|
|
|
|
|
|
# Commands
|
|
|
def cmd_fields(*args, **kwargs):
|
|
|
try:
|
|
|
path = os.path.join(BASE_PATH, args[0])
|
|
|
except IndexError:
|
|
|
path = '.'
|
|
|
|
|
|
ofile = os.path.join(path, 'field.dat')
|
|
|
|
|
|
fp = open(ofile, "w")
|
|
|
|
|
|
pattern = kwargs.get('pattern', 'C???')
|
|
|
subs = kwargs.get('subs', '')
|
|
|
last = kwargs.get('last', False)
|
|
|
|
|
|
get_grids = all_grids
|
|
|
|
|
|
if last:
|
|
|
get_grids = last_grid
|
|
|
|
|
|
|
|
|
for grid in get_grids(pattern=pattern, path=path):
|
|
|
print path, grid + subs
|
|
|
lz, m = load_axis('ez', grid + subs, path=path)
|
|
|
z, field_max, field_min = max_in_axis(lz, m, -1)
|
|
|
if grid != 'C000':
|
|
|
fp.write("%d\t%f\t%.18e\t%.18e\n" % (int(grid[1:]), z,
|
|
|
field_max, field_min))
|
|
|
|
|
|
print int(grid[1:]), z, field_max, field_min
|
|
|
fp.flush()
|
|
|
|
|
|
fp.close()
|
|
|
|
|
|
def cmd_fields_last(*args, **kwargs):
|
|
|
_d = kwargs.copy()
|
|
|
_d['last'] = True
|
|
|
|
|
|
cmd_fields(*args, **_d)
|
|
|
|
|
|
def cmd_posneg(*args, **kwargs):
|
|
|
try:
|
|
|
path = os.path.join(BASE_PATH, args[0])
|
|
|
except IndexError:
|
|
|
path = '.'
|
|
|
|
|
|
ofile = os.path.join(path, 'posneg.dat')
|
|
|
|
|
|
fp = open(ofile, "w")
|
|
|
pattern = kwargs.get('pattern', 'C???')
|
|
|
for grid in all_grids(pattern=pattern, path=path):
|
|
|
print path, grid
|
|
|
try:
|
|
|
neg_z, pos_z, neg_max, pos_max = posneg_fronts(grid, path=path)
|
|
|
|
|
|
if grid != 'C000':
|
|
|
fp.write("%d\t%f\t%f\t%f\t%f\n"
|
|
|
% (int(grid[1:]), neg_z, pos_z, neg_max, pos_max))
|
|
|
|
|
|
print int(grid[1:]), neg_z, pos_z, neg_max, pos_max
|
|
|
fp.flush()
|
|
|
except IOError:
|
|
|
print >> sys.stderr, "IOError exception in %s" % grid
|
|
|
|
|
|
fp.close()
|
|
|
|
|
|
|
|
|
def cmd_complete(*args, **kwargs):
|
|
|
try:
|
|
|
path = os.path.join(BASE_PATH, args[0])
|
|
|
except IndexError:
|
|
|
path = '.'
|
|
|
|
|
|
ofile = os.path.join(path, kwargs.get('ofile', 'analysis.dat'))
|
|
|
|
|
|
fp = open(ofile, "w")
|
|
|
pattern = kwargs.get('pattern', 'C???')
|
|
|
sub = kwargs.get('sub', '')
|
|
|
zmax = kwargs.get('zmax', None)
|
|
|
dt = kwargs.get('dt', 1.0)
|
|
|
|
|
|
try:
|
|
|
set_def_rescaling(kwargs['rescaling'])
|
|
|
except KeyError:
|
|
|
pass
|
|
|
|
|
|
|
|
|
|
|
|
fmt_str = None
|
|
|
|
|
|
varnames = ['t',
|
|
|
'Zc1', 'Zc2', 'Qmin', 'Qmax',
|
|
|
'Ze1', 'Ze2', 'Emin', 'Emax',
|
|
|
'V1', 'V2', 'R1', 'R2',
|
|
|
'Q', 'Q1', 'Q2']
|
|
|
|
|
|
fp.write("\t".join(varnames) + "\n")
|
|
|
|
|
|
old_c_fronts = None
|
|
|
pos_z0, neg_z0 = None, None
|
|
|
|
|
|
for grid in all_grids(pattern=pattern, path=path):
|
|
|
sgrid = grid + sub
|
|
|
clear_load_var_memoize()
|
|
|
print path, sgrid
|
|
|
try:
|
|
|
it = int(grid[1:]) * dt
|
|
|
itms = [it]
|
|
|
c_fronts = posneg_fronts(sgrid, path=path,
|
|
|
var='charge', zmax=zmax)
|
|
|
|
|
|
if pos_z0 is None or neg_z0 is None:
|
|
|
neg_z0, pos_z0 = [c_fronts[i] for i in 0, 1]
|
|
|
|
|
|
c_fronts[0] -= neg_z0
|
|
|
c_fronts[1] -= pos_z0
|
|
|
|
|
|
if old_c_fronts:
|
|
|
vminus, vplus = [(c_fronts[i] - old_c_fronts[i]) / (it - old_it)
|
|
|
for i in 0, 1]
|
|
|
else:
|
|
|
vminus, vplus = num.NAN, num.NAN
|
|
|
|
|
|
old_it = it
|
|
|
old_c_fronts = c_fronts
|
|
|
|
|
|
neg_r, pos_r = [fit_radius('charge', sgrid, s,
|
|
|
path=path, zmax=zmax)
|
|
|
for s in -1, 1]
|
|
|
|
|
|
#neg_diam, pos_diam = find_diameters('electrons', grid,
|
|
|
# c_fronts[0], c_fronts[1],
|
|
|
# path=path)
|
|
|
|
|
|
itms.extend(c_fronts)
|
|
|
|
|
|
e_fronts = posneg_fronts(sgrid, path=path, var='ez', zmax=zmax)
|
|
|
e_fronts[0] -= neg_z0
|
|
|
e_fronts[1] -= pos_z0
|
|
|
|
|
|
itms.extend(e_fronts)
|
|
|
itms.extend((vminus, vplus, neg_r, pos_r))
|
|
|
itms.append(integrate(sgrid, path=path, var='charge'))
|
|
|
|
|
|
# This was the old way of calculating the integrated charge:
|
|
|
# neg_zrange=(neg_z0 - 2*neg_r, neg_z0 + neg_r)
|
|
|
# pos_zrange=(pos_z0 - 2*pos_r, pos_z0 + pos_r)
|
|
|
#
|
|
|
# itms.append(integrate(sgrid, path=path, var='charge',
|
|
|
# sign=-1.0, threshold=0.5, zrange=neg_zrange))
|
|
|
# itms.append(integrate(sgrid, path=path, var='charge',
|
|
|
# sign=1.0, threshold=0.5, zrange=pos_zrange))
|
|
|
|
|
|
neg_zrange=(neg_z0 - neg_r, neg_z0 + 0.5 * neg_r)
|
|
|
pos_zrange=(pos_z0 - pos_r, pos_z0 + 0.5 * pos_r)
|
|
|
|
|
|
itms.append(integrate(sgrid, path=path, var='charge',
|
|
|
zrange=neg_zrange))
|
|
|
itms.append(integrate(sgrid, path=path, var='charge',
|
|
|
zrange=pos_zrange))
|
|
|
|
|
|
if fmt_str == None:
|
|
|
fmt_str = "\t".join("%g" for x in itms) + "\n"
|
|
|
|
|
|
if int(grid[1:4]) >= 2:
|
|
|
fp.write(fmt_str % tuple(itms))
|
|
|
fp.flush()
|
|
|
|
|
|
print "\t".join(str(x) for x in itms)
|
|
|
except IOError:
|
|
|
print >> sys.stderr, "IOError exception in %s" % grid
|
|
|
|
|
|
fp.close()
|
|
|
|
|
|
|
|
|
|
|
|
def cmd_diameters(runid, var, **kwargs):
|
|
|
try:
|
|
|
path = os.path.join(BASE_PATH, runid)
|
|
|
except IndexError:
|
|
|
path = '.'
|
|
|
|
|
|
ofile = os.path.join(path, 'diameters.dat')
|
|
|
fp = open(ofile, "w")
|
|
|
|
|
|
pattern = kwargs.get('pattern', 'C???')
|
|
|
for grid in all_grids(pattern=pattern, path=path):
|
|
|
try:
|
|
|
|
|
|
neg_z, pos_z, neg_max, pos_max = posneg_fronts(grid, path=path)
|
|
|
neg_diam, pos_diam = find_diameters(var, grid, neg_z, pos_z, path=path)
|
|
|
print grid
|
|
|
|
|
|
if grid != 'C000':
|
|
|
fp.write("%d\t%f\t%f\t%f\t%f\t%f\t%f\n"
|
|
|
% (int(grid[1:]), neg_z, pos_z,
|
|
|
neg_max, pos_max, neg_diam, pos_diam))
|
|
|
|
|
|
print int(grid[1:]), neg_z, pos_z, neg_max, pos_max, neg_diam,\
|
|
|
pos_diam
|
|
|
fp.flush()
|
|
|
except IOError:
|
|
|
print >> sys.stderr, "IOError exception in %s" % grid
|
|
|
|
|
|
fp.close()
|
|
|
|
|
|
|
|
|
def cmd_voltage(*args, **kwargs):
|
|
|
try:
|
|
|
path = os.path.join(BASE_PATH, args[0])
|
|
|
except IndexError:
|
|
|
path = '.'
|
|
|
|
|
|
ofile = os.path.join(path, 'voltage.dat')
|
|
|
|
|
|
fp = open(ofile, "w")
|
|
|
pattern = kwargs.get('pattern', 'C???')
|
|
|
for grid in all_grids(pattern=pattern, path=path):
|
|
|
try:
|
|
|
lz, m = load_axis('ez', grid, path=path)
|
|
|
voltage = - sum(m) * (lz[1] - lz[0])
|
|
|
s = "%d\t%f\n" % (int(grid[1:]), voltage)
|
|
|
fp.write(s)
|
|
|
fp.flush()
|
|
|
print s
|
|
|
|
|
|
except IOError:
|
|
|
print >> sys.stderr, "IOError exception in %s" % grid
|
|
|
|
|
|
|
|
|
def cmd_front(grid, ofile):
|
|
|
fp = open(ofile, "w")
|
|
|
|
|
|
id, grid = grid.split(':')
|
|
|
path = os.path.join(BASE_PATH, id)
|
|
|
|
|
|
print grid, path
|
|
|
|
|
|
front = find_front('charge', grid, path=path, fmin=0.0001, rmin=2.1, rmax=300,
|
|
|
zmin=0, zmax=2000, shift_origin=False,
|
|
|
pre_margin=200, post_margin=20)
|
|
|
for z, r, f in front:
|
|
|
fp.write("%f\t%f\t%f\n" % (z, r, f))
|
|
|
|
|
|
fp.close()
|
|
|
|
|
|
def cmd_fit_saffman(grid):
|
|
|
id, grid = grid.split(':')
|
|
|
return fit_saffman(id, grid)
|
|
|
|
|
|
def cmd_envelope(grid, ofile):
|
|
|
id, grid = grid.split(':')
|
|
|
path = os.path.join(BASE_PATH, id)
|
|
|
|
|
|
front = find_envelope(grid, path=path)
|
|
|
|
|
|
fp = open(ofile, "w")
|
|
|
for r, z in front:
|
|
|
fp.write("%f\t%f\n" % (r, z))
|
|
|
fp.close()
|
|
|
|
|
|
return front
|
|
|
|
|
|
def cmd_simple_front(var, grid):
|
|
|
id, grid = grid.split(':')
|
|
|
path = os.path.join(BASE_PATH, id)
|
|
|
|
|
|
front = simple_find_front(var, grid, path=path)
|
|
|
for r, z in front:
|
|
|
print r, z
|
|
|
|
|
|
return front
|
|
|
|
|
|
def cmd_iterfront(grid, n, ofile):
|
|
|
id, grid = grid.split(':')
|
|
|
path = os.path.join(BASE_PATH, id)
|
|
|
|
|
|
fp = open(ofile, "w")
|
|
|
n = int(n)
|
|
|
|
|
|
front = find_front_iter('eabs', grid, n, path=path,
|
|
|
fmin=0.001, rmin=1, rmax=300,
|
|
|
zmin=300, zmax=600, shift_origin=False,
|
|
|
pre_margin=250, post_margin=20)
|
|
|
for z, r, f in front:
|
|
|
fp.write("%f\t%f\t%f\n" % (z, r, f))
|
|
|
|
|
|
fp.close()
|
|
|
|
|
|
|
|
|
def cmd_max_off_axis(runid, var, sgn, **kwargs):
|
|
|
try:
|
|
|
path = os.path.join(BASE_PATH, runid)
|
|
|
except IndexError:
|
|
|
path = '.'
|
|
|
|
|
|
ofile = os.path.join(path, 'max-off-%s.dat' % var)
|
|
|
fp = open(ofile, "w")
|
|
|
|
|
|
pattern = kwargs.get('pattern', 'C???')
|
|
|
sgn = int(sgn)
|
|
|
for grid in all_grids(pattern=pattern, path=path):
|
|
|
try:
|
|
|
lr, lz, m = load_var(var, grid, path=path)
|
|
|
except IOError:
|
|
|
print >> sys.stderr, "IOError exception in %s" % grid
|
|
|
continue
|
|
|
|
|
|
r, z = max_off_axis(lr, lz, m, sgn)
|
|
|
if grid != 'C000':
|
|
|
fp.write("%d\t%f\t%f\n" % (int(grid[1:]), r, z))
|
|
|
|
|
|
print int(grid[1:]), r, z
|
|
|
fp.flush()
|
|
|
|
|
|
fp.close()
|
|
|
|
|
|
def cmd_shell():
|
|
|
from IPython.Shell import IPShellEmbed
|
|
|
ipshell = IPShellEmbed([])
|
|
|
ipshell()
|
|
|
|
|
|
def check_eval_float(option, opt, value):
|
|
|
eval_value = eval(value)
|
|
|
try:
|
|
|
return float(eval_value)
|
|
|
except ValueError:
|
|
|
raise OptionValueError(
|
|
|
"option %s: invalid float value: '%s' -> %r"
|
|
|
% (opt, value, eval_value))
|
|
|
|
|
|
|
|
|
|
|
|
class ExtOption(Option):
|
|
|
TYPES = Option.TYPES + ("eval_float",)
|
|
|
TYPE_CHECKER = Option.TYPE_CHECKER.copy()
|
|
|
TYPE_CHECKER["eval_float"] = check_eval_float
|
|
|
|
|
|
|
|
|
|
|
|
def main():
|
|
|
parser = OptionParser("", option_class=ExtOption)
|
|
|
|
|
|
parser.add_option('--dt', '', action='store',
|
|
|
dest='dt',
|
|
|
type='eval_float',
|
|
|
default=1.0,
|
|
|
help="dt")
|
|
|
|
|
|
parser.add_option('--with-units', '', action='store',
|
|
|
dest='units',
|
|
|
default=None,
|
|
|
help="Normalize with units")
|
|
|
|
|
|
|
|
|
parser.add_option('--pattern', '', action='store',
|
|
|
dest='pattern', type='str',
|
|
|
default='C???',
|
|
|
help='Pattern of the analyzed grids')
|
|
|
|
|
|
|
|
|
parser.add_option('--sub', '', action='store',
|
|
|
dest='sub', type='str',
|
|
|
default="",
|
|
|
help='Sub-grid postfix')
|
|
|
|
|
|
parser.add_option('--ofile', '', action='store',
|
|
|
dest='sub', type='str',
|
|
|
help='Output file')
|
|
|
|
|
|
|
|
|
parser.add_option('--zmax', '', action='store',
|
|
|
dest='zmax', type='eval_float',
|
|
|
default=None,
|
|
|
help='Sub-grid postfix')
|
|
|
|
|
|
parser.add_option('--rescaling', '', action='store',
|
|
|
dest='rescaling', type='str',
|
|
|
default=None,
|
|
|
help='Rescaling to use')
|
|
|
|
|
|
|
|
|
main = sys.modules[__name__]
|
|
|
|
|
|
opts, args = parser.parse_args()
|
|
|
|
|
|
try:
|
|
|
cmd = sys.argv[1].replace('-', '_')
|
|
|
process_method = getattr(main, 'cmd_%s' % cmd)
|
|
|
except AttributeError:
|
|
|
print "Dont know how to process command `%s`" % args[0]
|
|
|
sys.exit(-1)
|
|
|
|
|
|
process_method(*args[1:], **opts.__dict__)
|
|
|
|
|
|
|
|
|
|
|
|
if __name__ == "__main__":
|
|
|
main()
|