require 'narray'
class NArray
def _dump *ignored
Marshal.dump :typecode => typecode, :shape => shape, :data => to_s
end
def self._load buf
h = Marshal.load buf
typecode = h[:typecode]
shape = h[:shape]
data = h[:data]
to_na data, typecode, *shape
end
def inspect
#ep "called inspect"
"#{self.class}.to_narray(#{self.to_a.inspect})"
end
end
class GSL::Tensor
class << self
def method_missing(meth, *args)
#ep 'calling... ', meth, args
ans = new(NArray.send(meth, *args.reverse))
#ep 'got', ans
ans
end
end
attr_reader :narray
def self.alloc(*args)
new(NArray.float(*args.reverse))
end
def initialize(narray)
@narray = narray
end
def inspect
"GSL::Tensor.new(#{@narray.inspect})"
end
def [](*args)
#if args.inject(true){|b,i| b and i.kind_of? Integer}
#@narray[*args.reverse]
#else
#self.class.new(@narray[*args.reverse])
#end
case ans = @narray[*args.reverse]
when Numeric
ans
else
self.class.new(@narray[*args.reverse])
end
end
def []=(*args, value)
#def []=(*args)
#ep 'args', args, value
@narray[*args.reverse] = value
end
def shape
@narray.shape.reverse
end
def reshape!(*args)
#ep 'rags', args
@narray.reshape!(*args.reverse)
end
def to_a
@narray.transpose(*(0...@narray.shape.size).to_a.reverse).to_a
end
def transpose(*args)
self.class.new(@narray.transpose(*args))
end
def method_missing(meth, *args)
result = @narray.send(meth, *args.reverse)
if result.kind_of? NArray
self.class.new(result)
else
result
end
rescue NoMethodError
self.class.new(NMath.send(meth, @narray))
end
def iterate(&block)
shp = shape
cumul = 1
cumulshp = []
for i in 1..shape.size
cumulshp[shp.size-i] = cumul
cumul *= shp[shp.size-i]
end
#= shape.reverse.map{|dim| cumul*=(dim); cumul.to_i}.reverse
#ep cumulshp; gets
(cumulshp[0]*shp[0]).times do |n|
#indexes = cumulshp.reverse.map{|cumul| rem = n%cumul; n-=rem; rem}.reverse
indexes = cumulshp.map{|cumul| idx = (n/cumul).floor; n -= idx*cumul; idx}
yield(*indexes)
end
end
def iterate_row_maj(&block)
shp = shape
cumul = 1
cumulshp = []
for i in 0...shape.size
cumulshp[i] = cumul
cumul *= shp[i]
end
#= shape.reverse.map{|dim| cumul*=(dim); cumul.to_i}.reverse
#ep cumulshp; gets
(cumulshp[-1]*shp[-1]).times do |n|
#indexes = cumulshp.reverse.map{|cumul| rem = n%cumul; n-=rem; rem}.reverse
indexes = cumulshp.reverse.map{|cumul| idx = (n/cumul).floor; n -= idx*cumul; idx}.reverse
yield(*indexes)
end
end
end
class GSL::TensorComplex < GSL::Tensor
attr_reader :narray
def self.alloc(*args)
new(NArray.complex(*args.reverse))
end
def real
GSL::Tensor.new(@narray.real)
end
def abs
GSL::Tensor.new(@narray.abs)
end
end
class CodeRunner::Gs2
def gsl_tensor(name, options)
tensor = send((name.to_s+"_gsl_tensor").to_sym , options)
end
module GSLTensors
def moment_gsl_tensor(options)
if options[:t_index]
raise ArgumentError.new("Moments are not written out as a function of time currently")
#ep options; gets
raise CRFatal.new("write_phi_over_time is not enabled so this function won't work") unless @write_phi_over_time
arr = GSL::Tensor.new(netcdf_file.var(options[:field_name].to_s + '_t').get({'start' => [0,(options[:thetamin]||0),0,0, options[:t_index] - 1], 'end' => [-1,(options[:thetamax]||-1),(options[:nakx]||0)-1,(options[:naky]||0)-1, options[:t_index] - 1]}))
#ep 'arr.shape', arr.shape
arr.reshape!(*arr.shape.slice(1...arr.shape.size))
else
arr = GSL::Tensor.new(netcdf_file.var(options[:moment_name]).get({'start' => [0,(options[:thetamin]||0),0,0,options[:species_element]], 'end' => [-1,(options[:thetamax]||-1),(options[:nakx]||0)-1,(options[:naky]||0)-1,options[:species_element]]}))
#ep 'arr.shape', arr.shape
end
arr.reshape!(*arr.shape.slice(1...arr.shape.size))
arr[0, true, true, true] = 0.0 if options[:no_zonal]
#arr = arr[options[:nakx] ? 0...options[:nakx] : true, options[:naky] ? 0...options[:naky] : true, true, true] if options[:nakx] or options[:naky]
return arr
end
def field_netcdf_name(field_name, time_varying = false)
#p field_name.to_s
name = case field_name.to_s
when /phi/
time_varying ? 'phi_t' : 'phi'
when /density/
time_varying ? 'ntot_t' : 'density'
when /apar/
time_varying ? 'apar_t' : 'apar'
else
raise "Unknown field name: #{field_name}"
end
#p name
return name
end
def field_species_element(options)
case options[:field_name].to_s
when /density/
options.convert_to_index(self, :species)
#ep 'options', options
options[:species_index] - 1
else
nil
end
end
def field_gsl_tensor(options)
species_element = field_species_element(options)
#ep 'species_element', species_element
if options[:t_index]
#ep options; gets
#raise CRFatal.new("write_phi_over_time is not enabled so this function won't work") unless @write_phi_over_time
arr = GSL::Tensor.new(netcdf_file.var(field_netcdf_name(options[:field_name], true)).get({'start' => [0,(options[:thetamin]||0),0,0, species_element, options[:t_index] - 1].compact, 'end' => [-1,(options[:thetamax]||-1),(options[:nakx]||0)-1,(options[:naky]||0)-1, species_element, options[:t_index] - 1].compact}))
#ep 'arr.shape', arr.shape
arr.reshape!(*arr.shape.slice(1...arr.shape.size))
else
arr = GSL::Tensor.new(netcdf_file.var(field_netcdf_name(options[:field_name])).get({'start' => [0,(options[:thetamin]||0),0,0, species_element].compact, 'end' => [-1,(options[:thetamax]||-1),(options[:nakx]||0)-1,(options[:naky]||0)-1, species_element].compact}))
#ep 'arr.shape', arr.shape
end
if species_element
arr.reshape!(*arr.shape.slice(1...arr.shape.size))
end
if options[:interpolate_x]
shape = arr.narray.shape
#p 'shape', shape
shape[2] = (shape[2]-1)*options[:interpolate_x] + 1
#p shape
arr = GSL::Tensor.new(arr.narray.expand(*shape, 0.0))
end
if options[:interpolate_y]
shape = arr.narray.shape
#p 'shape', shape
shape[3] = (shape[3]-1)*options[:interpolate_y] + 1
#p shape
arr = GSL::Tensor.new(arr.narray.expand(*shape, 0.0))
end
if gryfx? and options[:periodic]
shape = arr.narray.shape
shape[1]+=1
arr = GSL::Tensor.new(arr.narray.expand(*shape, 0.0))
shpe = arr.shape
for i in 0...shpe[0]
for j in 0...shpe[1]
for r in 0...shpe[3]
arr[i, j, -1, r] = arr[i, j, 0, r]
end
end
end
end
arr[0, true, true, true] = 0.0 if options[:no_zonal]
#arr = arr[options[:nakx] ? 0...options[:nakx] : true, options[:naky] ? 0...options[:naky] : true, true, true] if options[:nakx] or options[:naky]
return arr
end
# Returns a rank 3 tensor which is the real potential (i.e. Fourier transformed from the GS2 output) as a function of the y index, the x index and the theta index.
def phi_real_space_gsl_tensor(options)
return field_real_space_gsl_tensor(options.absorb(field_name: :phi))
end
def field_real_space_gsl_tensor(options)
fieldc = field_gsl_tensor_complex(options)
shape = fieldc.shape
workspacex = GSL::Vector::Complex.alloc(shape[1])
workspacey = GSL::Vector.alloc(shape[0]*2-2+shape[0]%2)
field_real_space = GSL::Tensor.alloc(workspacey.size, shape[1], shape[2])
for j in 0...shape[2] #theta
for i in 0...shape[0] #ky
#narr = fieldc[i, true, j]
for k in 0...shape[1]
workspacex[k] = GSL::Complex.alloc(fieldc[i,k,j].real, fieldc[i,k,j].imag)
end
workspacex = workspacex.backward
for k in 0...shape[1]
fieldc[i,k,j] = Complex(*workspacex[k].to_a)
end
end
for k in 0...shape[1] #kx
m = 0
for i in 0...shape[0] #ky
workspacey[m] = fieldc[i,k,j].real
m+=1
next if i==0 or (shape[0]%2==0 and i == shape[0]/2 + 1)
workspacey[m] = fieldc[i,k,j].imag
m+=1
end
workspacey = workspacey.backward
for i in 0...workspacey.size
field_real_space[i,k,j] = workspacey[i]
end
end
end
shp = field_real_space.shape
#ep options
field_real_space = field_real_space[options[:ymin]||0..options[:ymax]||(shp[0]-1), options[:xmin]||0..options[:xmax]||(shp[1]-1), true]
if kint = options[:interpolate_theta]
shape = field_real_space.shape
new_shape = shape.dup
new_shape[-1] = ((shape[-1]-1)*kint+1)
field_real_space_new = GSL::Tensor.float(*new_shape)
#p shape,new_shape
for i in 0...(new_shape[0])
for j in 0...(new_shape[1])
field_real_space_new[i,j, new_shape[-1]-1] = field_real_space[i,j,shape[-1]-1] # set the endpoint
for k in 0...(new_shape[-1]-1)
km = k%kint
frac = km.to_f/kint.to_f
#kold = (k-km)/(new_shape[-1]-1)*(shape[-1]-1)
kold = (k-km)/kint
#ep ['k', k, 'kold', kold]
field_real_space_new[i,j, k] = field_real_space[i,j, kold] * (1.0-frac) + field_real_space[i,j, kold+1] * frac
end
end
end
field_real_space = field_real_space_new
end
return field_real_space
end
def field_real_space_gsl_tensor_2(options)
field = field_gsl_tensor(options)
field_narray = field.narray
shape = field.shape
workspacex = GSL::Vector::Complex.alloc(shape[1])
workspacey = GSL::Vector.alloc(shape[0]*2-2+shape[0]%2)
field_real_space = GSL::Tensor.alloc(workspacey.size, shape[1], shape[2])
field_real_space_narray = field_real_space.narray
for j in 0...shape[2] #theta
for i in 0...shape[0] #ky
#narr = fieldc[i, true, j]
for k in 0...shape[1]
workspacex[k] = GSL::Complex.alloc(field_narray[0,j,k,i], field_narray[1,j,k,i])
end
workspacex = workspacex.backward
for k in 0...shape[1]
field_narray[0,j,k,i] = workspacex[k].real
field_narray[1,j,k,i] = workspacex[k].imag
end
end
for k in 0...shape[1] #kx
m = 0
for i in 0...shape[0] #ky
workspacey[m] = field_narray[0,j,k,i]
m+=1
next if i==0 or (shape[0]%2==0 and i == shape[0]/2 + 1)
workspacey[m] = field_narray[1,j,k,i]
m+=1
end
workspacey = workspacey.backward
for i in 0...workspacey.size
field_real_space_narray[j,k,i] = workspacey[i]
end
end
end
shp = field_real_space.shape
#p 'test', field_real_space[0,2,3]
#ep options
field_real_space = field_real_space[options[:ymin]||0..options[:ymax]||(shp[0]-1), options[:xmin]||0..options[:xmax]||(shp[1]-1), true]
#p 'test2', field_real_space[0,2,3]
if kint = options[:interpolate_theta]
shape = field_real_space.shape
new_shape = shape.dup
new_shape[-1] = ((shape[-1]-1)*kint+1)
field_real_space_new = GSL::Tensor.float(*new_shape)
field_real_space_new_narray = field_real_space_new.narray
#p shape,new_shape
for i in 0...(new_shape[0])
for j in 0...(new_shape[1])
field_real_space_new_narray[new_shape[-1]-1, j, i] = field_real_space_narray[shape[-1]-1, j, i] # set the endpoint
for k in 0...(new_shape[-1]-1)
km = k%kint
frac = km.to_f._orig_div(kint.to_f)
#kold = (k-km)/(new_shape[-1]-1)*(shape[-1]-1)
kold = (k-km)._orig_div(kint)
#ep ['k', k, 'kold', kold]
field_real_space_new_narray[k,j,i] = field_real_space_narray[kold,j,i]._orig_mul(1.0-frac) + field_real_space_narray[kold+1,j,i]._orig_mul(frac)
#if (i==0 and j==2 and k==3)
#p ['frac', frac]
#end
end
end
end
field_real_space = field_real_space_new
end
#p field_real_space_new.shape;
return field_real_space
end
def apar_gsl_tensor(options)
return GSL::Tensor.new(netcdf_file.var('apar').get)
end
def bpar_gsl_tensor(options)
return GSL::Tensor.new(netcdf_file.var('bpar').get)
end
# Order is R0,Z0,a0,Rprim,Zprim,aprim
def geometric_factors_gsl_tensor(options)
#ops = options.dup; ops.delete :phi
#ep ops; gets
case @equilibrium_option
when "s-alpha"
return geometric_factors_salpha_gsl_tensor(options)
else
theta_vec = gsl_vector(:theta, options)
factors = GSL::Tensor.alloc(6,theta_vec.size)
values = File.read(options[:geometry_file]||"#@directory/#@run_name.g").split(/\s*\n\s*/)
3.times{values.shift}
values = values.map{|str| str.split(/\s+/).map{|s| s.to_f}}.transpose
#ep values
shape = factors.shape
for i in 0...shape[0]
unless options[:interpolate_theta]
for j in 0...shape[1]
factors[i,j] = values[i+1][j]
end
else
opts = options.dup
opts[:interpolate_theta] = nil
theta_vec_short = gsl_vector(:theta, {})
interp = GSL::ScatterInterp.alloc(:linear, [values[0], values[i+1].to_gslv], true)
for j in 0...theta_vec.size
factors[i,j] = interp.eval(theta_vec[j])
end
end
end
#ep factors
return factors
end
end
# Order is R0,Z0,a0,Rprim,Zprim,aprim
def geometric_factors_salpha_gsl_tensor(options)
raise "Please specify options[:Rgeo]" unless options[:Rgeo]
theta_vec = gsl_vector(:theta, options)
factors = GSL::Tensor.alloc(6,theta_vec.size)
q_actual = options[:q_actual]
pka = epsl/q_actual
#pka = @pk||2*@kp
for i in 0...theta_vec.size
theta = theta_vec[i]
c = Math.cos(theta)
s = Math.sin(theta)
factors[0,i] = options[:Rgeo]*(1.0 + eps * c + eps * (shift||0))
factors[1,i] = options[:Rgeo] * eps * s
factors[2,i] = -epsl/pka * theta - eps * epsl/pka * s
factors[3,i] = c * options[:Rgeo] * eps / 2
factors[4,i] = s * options[:Rgeo] * eps / 2
factors[5,i] = - theta * epsl**2 / 2 / pka / eps * (shat||0) - epsl**2 / 2 / pka * s
end
return factors
end
private :geometric_factors_salpha_gsl_tensor
# Returns a rank 2 tensor, which gives, as a function of the x index j and the theta index k, the y index nearest to a poloidal plane at angle options[:torphi] is the torus was filled with periodic copies of the flux surface. Used for making cross sections at a constant toroidal angle.
def constant_torphi_surface_gsl_tensor(options)
ops = options.dup
IRRELEVANT_INDICES.each{|v| ops.delete(v)}
return cache[[:constant_torphi_surface_gsl_tensor, ops]] if cache[[:constant_torphi_surface_gsl_tensor, ops]]
correct_3d_options(options)
if options[:torphi]
torphiout = options[:torphi]
else
raise 'Need to specify a torphi value'
end
cyls = cylindrical_coordinates_gsl_tensor(options.absorb({extra_points: :y}))
shpc = cyls.shape
factors = geometric_factors_gsl_tensor(options)
#ep shpc, 'shpc'
#xsize = case shpc[2]
yvec = gsl_vector('y', options)
#ep yvec.to_a ; gets
x = gsl_vector('x', options)
dy = yvec[1] - yvec[0]
torphi_const = GSL::Tensor.int(shpc[2], shpc[3]) # don't include extra x point
xfac = 1.0 / options[:rho_star_actual]
yfac = options[:rhoc_actual] / options[:q_actual] / options[:rho_star_actual]
#coordinates[2,i,j,k] = y[i] / yfac - factors[2,k] - x[j]/xfac*factors[5,k] # phi
twopi = Math::PI*2
for j in 0...shpc[2]
for k in 0...shpc[3]
y = yfac * (torphiout + factors[2,k] + x[j]/xfac*factors[5,k])
if options[:no_copies]
i = (y/dy).floor
else
i = (y/dy).floor % yvec.size
end
torphi_const[j,k] = i
end
end
return torphi_const
#ep torphi_const; gets
end
#def constant_torphi_surface_gsl_tensor2(options)
#ops = options.dup
#IRRELEVANT_INDICES.each{|v| ops.delete(v)}
#return cache[[:constant_torphi_surface_gsl_tensor, ops]] if cache[[:constant_torphi_surface_gsl_tensor, ops]]
#torphiout = options[:torphi]
#correct_3d_options(options)
#cyls = cylindrical_coordinates_gsl_tensor(options.absorb({extra_points: :y}))
#shpc = cyls.shape
##ep shpc, 'shpc'
##xsize = case shpc[2]
#torphi_const = GSL::Tensor.int(shpc[2], shpc[3]) # don't include extra x point
##ep torphi_const; gets
#y = gsl_vector('y', options)
#lastbracketed = nil
#lastj = -1
#for k in 0...shpc[3] # theta index
#for j in 0...(shpc[2]) # x index
#deltorphi = cyls[2,shpc[1]-1,j,k] - cyls[2,0,j,k]
#raise "Periodicity not satisfied: #{(2*Math::PI/deltorphi+1.0e-8)%1.0}, #{(2*Math::PI/deltorphi+1.0e-5)}" unless ((2.0*Math::PI/deltorphi)+1.0e-8)%1.0 < 1.0e-5
#m3 = (torphiout)%deltorphi
#for i in 0...shpc[1]-1 # y index, excluding periodic extra point
##p i
#torphi1 = cyls[2,i,j,k]
#torphi2 = cyls[2,i+1,j,k]
##ep cyls[2,true,j,k].to_a
#m1 = (torphi1 )%deltorphi
#m2 = (torphi2 )%deltorphi
#bracketed = ((m1-m3).abs < 1.0e-4) || (
#(m2-m3.abs) > 1.0e-4 &&
#(m2 - m3) *
#(m1 - m3) *
#(m2 - m1) *
#(torphi2 - torphi1) < 0)
##p 'n0', (2*Math::PI/deltorphi).round
#bracketed2 = (2*Math::PI/deltorphi + 1).round.times.inject(false) do |b,n|
#epsn = 1.0e-4
#eps2 = 1.0e-4
#upp = torphiout + deltorphi * n
#lwr = torphiout - deltorphi * n
##measure = ((torphi1 < upp or (torphi1 - upp).abs < epsn) and upp+epsn < torphi2) or ((torphi1 < lwr or (torphi1-lwr).abs < eps) and lwr + eps < torphi2)
#a1 = a2 = a3 = b1 = b2 = b3 = 'q'
##measure = ((
###a1=((torphi1-upp).abs < (torphi2-upp).abs) and
##(a1=((torphi2-upp).abs > 1.0e-7)) and
##(a2=((torphi1 < upp or (torphi1 - upp).abs < epsn))) and
##(a3=(upp < torphi2))
##) or (
###b1=((torphi1-lwr).abs < (torphi2-lwr).abs) and
##(b1=((torphi2-lwr).abs > 1.0e7)) and
##(b2=(torphi1 < lwr or (torphi1-lwr).abs < epsn)) and
##(b3 = (lwr < torphi2))
##))
#a1=((torphi2-upp).abs > eps2)
#a2=((torphi1 < upp or (torphi1 - upp).abs < epsn))
#a3=(upp < torphi2)
#b1=(torphi2-lwr).abs > eps2
#b2=(torphi1 < lwr or (torphi1-lwr).abs < epsn)
#b3 = (lwr < torphi2)
#measure = ((a1 and a2 and a3) or (b1 and b2 and b3))
##p 'measure', measure, [torphi1, torphi2, upp, lwr , i,j,k, y[i], y[(i+1)%y.size], deltorphi, n, a1, a2, a3, b1, b2, b3] if measure and j==0 #; gets if [j,k] == [5,8] # if measure
#b or measure
#end
##bracketed = bracketed2
#raise "Measures don't agree #{bracketed}, #{bracketed2}" unless bracketed2 == bracketed
##d2 = torphi2 - torphiout
##d1 = torphiout - torphi1
#if bracketed
#raise "Doubled up" if lastbracketed == [j,k]
#raise "Missed: #{j},#{k}, #{lastbracketed.inspect} #{[j-1,k].inspect} #{[shpc[2]-1, k-1].inspect}" unless lastbracketed == [j-1,k] or lastbracketed == [shpc[2]-1, k-1] if lastbracketed
#torphi_const[j,k] = i
#lastbracketed = [j,k]
#lastj = j
#end
#end # y loop
#end # x loop
#end # theta loop
##torphi_const2 = constant_torphi_surface_gsl_tensor2(options)
##p 'the same? ', torphi_const2.to_a, torphi_const.to_a, torphi_const2 == torphi_const
##exit
##exit
#cache[[:constant_torphi_surface_gsl_tensor, ops]] = torphi_const
## save the run to save the hard_cache
#return torphi_const
#end
FIELD_VALUES = [:phi, :density, :apar, :bpar]
TRIVIAL_INDICES = [:graphkit_name]
TIME_VARYING_INDICES = [:t_index, :begin_element, :end_element, :frame_index, :t_index_window]
IRRELEVANT_INDICES = FIELD_VALUES + TRIVIAL_INDICES + TIME_VARYING_INDICES
# Adjust n0, rho_star_actual and q_actual to ensure periodicity
#
def correct_3d_options(options)
raise "Please specify options[:rho_star] or options[:n0]" unless options[:rho_star] or options[:n0]
case @equilibrium_option
when "s-alpha"
qinp = epsl / (pk||2*kp)
#xfac = @epsl**4/options[:rho_star]/4/pka**2/@eps**2
#xfac_geo = 1
#yfac = 1/options[:rho_star]/@epsl*2*pka*@eps
#yfac_geo = 2*pka*@eps/@epsl**2
#yfac_geo = 2*pka*@eps/@epsl**2
options[:rhoc_actual] =rhoc = 2 * eps / epsl
else
options[:rhoc_actual] = rhoc = @rhoc
qinp = @qinp
end
#eputs "Checking that rho_star and q satisfy periodicity..."
rho_star_inp = options[:rho_star]
y = gsl_vector('y', options)
ly = (y[1]-y[0]) * (y.size)
n0_fac = 2.0*Math::PI * rhoc / ly
n0_inp = options[:n0] || n0_fac / qinp / rho_star_inp
if n0_inp%1.0==0.0
n0 = n0_inp
else
#eputs "Input n0 is equal to #{n0_inp}..."
n0 = n0_inp.ceil
#eputs "Set n0 to #{n0}..."
end
if (qinp*n0)%1.0==0.0
q_actual = qinp
else
q_actual = (qinp*n0).round.to_f/n0
#eputs "Set q to #{q_actual}..."
end
options[:q_actual] = q_actual
unless options[:rho_star_actual] and options[:rho_star_actual] == n0_fac/n0/q_actual
#eputs "Adjusting rho_star to satisfy periodicity ..."
options[:rho_star_actual] = n0_fac/n0/q_actual
#eputs "Set rhostar to #{options[:rho_star_actual]}..."
#eputs "Note... to avoid adjustment of q specify n0 as an input rather than rho_star. Make sure that n0 is an integer and n0 * q is an integer."
end
end
# Return a rank 4 tensor which give
# cylindrical coordinates R,Z,torphi as a function
# of gs2 coordinates y, x, theta.
#
# a = cylindrical_coordinates_gsl_tensor(options)
#
# # pseudocode
# R(y[i], x[j], theta[k]) = a[0,i,j,k]
# Z(y[i], x[j], theta[k]) = a[1,i,j,k]
# torphi(y[i], x[j], theta[k]) = a[2,i,j,k]
def cylindrical_coordinates_gsl_tensor(options)
ops = options.dup
(IRRELEVANT_INDICES + [:torphi, :torphi_values]).each{|v| ops.delete(v)}
return cache[[:cylindrical_coordinates_gsl_tensor, ops]] if cache[[:cylindrical_coordinates_gsl_tensor, ops]]
#ep ops; gets
#options = options.dup
x = gsl_vector('x', options)
y = gsl_vector('y', options)
ly = 2*Math::PI*y0#(y[1]-y[0]) * (y.size)
if [true,:x].include? options[:extra_points]
ep "Extending x..."
x = x.connect([2*x[-1] - x[-2]].to_gslv).dup
end
if [true,:y].include? options[:extra_points]
ep "Extending y..."
y = y.connect([2.0*y[-1] - y[-2]].to_gslv).dup
raise "ly corrected incorrectly #{ly},#{y[-1]},#{y[0]},#{y[-1]-y[0]}" unless (ly-(y[-1] - y[0])).abs / ly.abs < 1.0e-6
end
x = x.subvector(options[:xmin]||0, (options[:xmax]||x.size-1) - (options[:xmin]||0) + 1).dup # if options[:xout] and options[:xin]
y = y.subvector(options[:ymin]||0, (options[:ymax]||y.size-1) - (options[:ymin]||0) + 1).dup # if options[:yout] and options[:yin]
y = y + options[:ncopy] * ly if options[:ncopy]
theta = gsl_vector('theta', options)
correct_3d_options(options)
rhoc = options[:rhoc_actual]
q_actual = options[:q_actual]
xfac = 1.0 / options[:rho_star_actual]
yfac = rhoc / q_actual / options[:rho_star_actual]
factors = geometric_factors_gsl_tensor(options)
coordinates = GSL::Tensor.alloc(3, y.size, x.size, theta.size)
for i in 0...y.size
for j in 0...x.size
for k in 0...theta.size
coordinates[0,i,j,k] = factors[0,k] + x[j]/xfac*factors[3,k] # R
coordinates[1,i,j,k] = factors[1,k] + x[j]/xfac*factors[4,k] # Z
coordinates[2,i,j,k] = y[i] / yfac - factors[2,k] - x[j]/xfac*factors[5,k] # phi
if gs2f = options[:gs2_coordinate_factor]
rgs2 = (x[j]**2 + y[i]**2 )**0.5*(1.0 + 2.0 * Float::EPSILON)
if rgs2 < 1.0e-8
phigs2 = 0
else
phigs2 = Math.acos(x[j]/rgs2)
end
coordinates[0,i,j,k] = rgs2 * gs2f + coordinates[0,i,j,k] * (1.0-gs2f)
coordinates[1,i,j,k] = theta[k] * gs2f + coordinates[1,i,j,k] * (1.0-gs2f)
coordinates[2,i,j,k] = phigs2 * gs2f + coordinates[2,i,j,k] * (1.0-gs2f)
end
end
end
end
case tp = options[:toroidal_projection]
when Numeric
coordinates[2, false] = tp
end
cache[[:cylindrical_coordinates_gsl_tensor, ops]] = coordinates
#save # save the run to save the hard_cache
return coordinates
end
# Return a rank 4 tensor which give
# cartesian coordinates X,Y,Z as a function
# of gs2 coordinates y, x, theta.
#
# a = cartesian_coordinates_gsl_tensor(options)
#
# # pseudocode
# X(y[i], x[j], theta[k]) = a[0,i,j,k]
# Y(y[i], x[j], theta[k]) = a[1,i,j,k]
# Z(y[i], x[j], theta[k]) = a[2,i,j,k]
def cartesian_coordinates_gsl_tensor(options)
cyl = cylindrical_coordinates_gsl_tensor(options)
shape = cyl.shape
cart = GSL::Tensor.alloc(*shape)
for i in 0...shape[1]
for j in 0...shape[2]
for k in 0...shape[3]
r = cyl[0,i,j,k]
z = cyl[1,i,j,k]
phi = cyl[2,i,j,k]
#cart[0,i,j,k] = r # Y
cart[0,i,j,k] = r*Math.cos(phi) # X
#cart[1,i,j,k] = phi # X
cart[1,i,j,k] = r*Math.sin(phi) # Y
cart[2,i,j,k] = z
end
end
end
cart
end
end #module
include GSLTensors
module GSLComplexTensors
def phi_gsl_tensor_complex(options)
return field_gsl_tensor_complex(options.absorb({field_name: :phi}))
end
def field_gsl_tensor_complex(options)
field = field_gsl_tensor(options)
fieldc = GSL::TensorComplex.alloc(*field.shape.slice(0..2))
nac = fieldc.narray
na = field.narray
for i in 0...field.shape[0]
for j in 0...field.shape[1]
for k in 0...field.shape[2]
nac[k,j,i] = Complex(na[0,k,j,i],na[1,k,j,i])
end
end
end
return fieldc
end
end #module
include GSLComplexTensors
end #class