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