"""
infinite_fourier.py
====================================
Module for calculating :math:`C_I` using linear algebra,
based on components generated by the method of recursive images
using fourier series to obtain :math:`V_{I\infty}`
"""
import scipy
import scipy.special
import numpy as np
import numpy.linalg
eps0= 8.854187817*10**-12

class layer:
	def __init__(self,eps_x,eps_y,t):
		"""
		layer class

		Parameters
		----------
		eps_x : float
			relative dielectric constant in the x-direction
		eps_ys : layer object
			relative dielectric constant in the y-direction
		t : float
			thickness of layer
		max_fourier_n : int
			max number of fourier compoents in model
		"""
		self.eps_x=eps_x
		self.eps_y=eps_y
		if self.eps_y>0:
			self.eps_r=(self.eps_x/self.eps_y)**0.5
		else:
			self.eps_r=1
		self.t=t
		self.t_eff=t*self.eps_r
		self.interface_pos_dir=None
		self.interface_neg_dir=None
class interface:
	"""
	Interface class

	Parameters
	----------
	layer_pos_dir : layer object
		The layer in positive y-direction
	layer_neg_dir : layer object
		The layer in negative y-direction
	y : float
		y-coordinate of the layer
	max_fourier_n : int
		max number of fourier compoents in model
	"""

	def __init__(self,layer_pos_dir,layer_neg_dir,y,max_fourier_n):
		self.layer_pos_dir=layer_pos_dir
		self.layer_neg_dir=layer_neg_dir
		# calculate regular (not effective) r and t
		eps1=np.sqrt(layer_pos_dir.eps_y*layer_pos_dir.eps_x)
		eps2=np.sqrt(layer_neg_dir.eps_y*layer_neg_dir.eps_x)
		self.r_pos_dir=(eps2-eps1)/(eps1+eps2)
		self.t_pos_dir=self.r_pos_dir+1
		self.r_neg_dir=-self.r_pos_dir
		self.t_neg_dir=self.r_neg_dir+1
		self.y=y
		# set all effective coefficients to None, are calculated when the get functions are called
		self.r_pos_dir_eff=np.full((max_fourier_n), None)
		self.t_pos_dir_eff=np.full((max_fourier_n), None)
		self.r_neg_dir_eff=np.full((max_fourier_n), None)
		self.t_neg_dir_eff=np.full((max_fourier_n), None)
	def get_r_pos_dir_eff(self,n):
		"""
		Calculate :math:`\\r_{eff}` in the positive direction

		Parameters
		----------
		n : int
			Forier compoent
		"""
		if self.r_pos_dir_eff[n]==None:
			if self.layer_pos_dir.interface_pos_dir==None:
				# if this is the last material, it has infinite extension, and it is easy
				self.r_pos_dir_eff[n]=self.r_pos_dir
				self.t_pos_dir_eff[n]=self.t_pos_dir
			else:
				# if this is not the last material, calculate effective coefficients r and t
				rup=self.layer_pos_dir.interface_pos_dir.get_r_pos_dir_eff(n)
				t=self.layer_pos_dir.t_eff
				exp=np.exp(-np.pi*2*t*(2*n+1))
				factor=1/(1-self.r_neg_dir*rup*exp) #== sum {m=0->inf} (self.r_neg_dir*rup*exp)^m
				self.r_pos_dir_eff[n]=self.r_pos_dir+self.t_pos_dir*rup*self.t_neg_dir*exp*factor
				self.t_pos_dir_eff[n]=self.t_pos_dir+self.t_pos_dir*rup*self.r_neg_dir*exp*factor
		return self.r_pos_dir_eff[n]
	def get_t_pos_dir_eff(self,n):
		"""
		Calculate :math:`\\t_{eff}` in the positive direction

		Parameters
		----------
		n : int
			Forier compoent
		"""
		if self.t_pos_dir_eff[n]==None:
			self.get_r_pos_dir_eff(n)
		return self.t_pos_dir_eff[n]
	def get_r_neg_dir_eff(self,n):
		"""
		Calculate :math:`\\r_{eff}` in the negative direction

		Parameters
		----------
		n : int
			Forier compoent
		"""
		if self.r_neg_dir_eff[n]==None:
			if self.layer_neg_dir.interface_neg_dir==None:
				# if this is the first material, it has infinite extension, and it is easy
				self.r_neg_dir_eff[n]=self.r_neg_dir
				self.t_neg_dir_eff[n]=self.t_neg_dir
			else:
				# if this is not the last material, calculate effective coefficients r and t
				rne=self.layer_neg_dir.interface_neg_dir.get_r_neg_dir_eff(n)
				t=self.layer_neg_dir.t_eff
				exp=np.exp(-np.pi*2*t*(2*n+1))
				factor=1/(1-self.r_pos_dir*rne*exp)
				self.r_neg_dir_eff[n]=self.r_neg_dir+self.t_neg_dir*rne*self.t_pos_dir*exp*factor
				self.t_neg_dir_eff[n]=self.t_neg_dir+self.t_neg_dir*rne*self.r_pos_dir*exp*factor
		return self.r_neg_dir_eff[n]
	def get_t_neg_dir_eff(self,n):
		"""
		Calculate :math:`\\t_{eff}` in the negative direction

		Parameters
		----------
		n : int
			Forier compoent
		"""
		if self.t_neg_dir_eff[n]==None:
			self.get_r_neg_dir_eff(n)
		return self.t_neg_dir_eff[n]

class single_recursive_images:
	"""
		A class that houses the potential, capacitance and electric fields from the method of recursive images.
		The values obtained are unphysical, as the potential at the electrodes is not constant.

		Parameters
		----------
		eta : float
			cover fraction of the electrodes
		interface_of_electrodes : int
			interface for the electrodes
		thickness_of_layers : list of floats
			thicknesses of the layers, this list will be 2 shorter than eps_x_of_layers and eps_y_of_layers, as the outermost layers have no defined thickness, but are infinite
		eps_x_of_layers : list of floats
			in-plane dielectric constant of the layers
		eps_y_of_layers : list of floats
			out-of-plane dielectric constant of the layers
		max_fourier_n : int, optional
			maximum number of fourier components to considder, defaults to 120
		accuracy_limit : float, optional
			reflections with  less than accuracy_limit are ignored, defaults to 10**-15
		inherit: single_recursive_images object or None, optional
			If not None, will inherit layers and interfaces from the pobject indicated.
			If None, new layer and interface objects will be created. Defaults to None
		"""
	_A_ditcts=dict()#adict[k][n]=A_{2n+1}
	_V_dicts=dict()#Vdict[k][x]=V_{I,inf}(x,0)
	#Common variables
	#----------------
	#_A_ditcts : dictionary
	#	:math:`A(\\eta)_{2n+1}` is stored here so it does not need to be recalculated for later iterations
	#_V_dicts : dictionary
	#	:math:`V(\\eta)_{2n+1}` is stored here so it does not need to be recalculated for later iterations
	def __init__(self,eta,interface_of_electrodes,thickness_of_layers,eps_x_of_layers,eps_y_of_layers,max_fourier_n=120,accuracy_limit=10**-15, inherit=None, hybrid=True):
		"""
		Returns
		-------
		class object
			"""

		self.hybrid=hybrid
		# declare the layers and interfaces, if not inherited
		if inherit==None:
			self.layers=[]
			for e11,e33,T in zip(eps_x_of_layers,eps_y_of_layers,[np.inf]+list(thickness_of_layers)+[np.inf]):
				self.layers.append(layer(e11,e33,T))
			self.interfaces=[]
			for i, _ in enumerate(self.layers[0:-1]):
				if i==0:
					z=0
				else:
					z=sum(thickness_of_layers[0:i])
				self.interfaces.append(interface(self.layers[i+1],self.layers[i],z,max_fourier_n))
				self.layers[i].interface_pos_dir=self.interfaces[-1]
				self.layers[i+1].interface_neg_dir=self.interfaces[-1]
		else:
			self.layers=inherit.layers
			self.interfaces=inherit.interfaces
		self.eta=np.array(eta)
		self.interface_of_electrodes=interface_of_electrodes
		self.number_of_layers=len(self.layers)
		self.max_fourier_n=max_fourier_n
		self.accuracy_limit=accuracy_limit
		self.k=np.sin(np.pi/2*self.eta)
		self.Kk= scipy.special.ellipk(float(self.k**2))
		self.Kpk= scipy.special.ellipk(1-float(self.k**2))
		if not self.k in self._A_ditcts:
			self._A_ditcts[self.k]=dict()
		self.adict=self._A_ditcts[self.k] #dict()#adict[n]=A_{2n+1}
		if self.hybrid:
			if not self.k in self._V_dicts:
				self._V_dicts[self.k]=dict()
			self.Vdict=self._V_dicts[self.k] #dict()#adict[n]=A_{2n+1}
		'''if not self.k in self._V_dicts:
			self._V_dicts[self.k]=dict()
		self.Vdict=self._V_dicts[self.k] #dict()#adict[n]=A_{2n+1}'''
		self.tree=[]
		self.P_field_pos=np.full((max_fourier_n), None)#r_pos_dir)
		self.P_field_neg=np.full((max_fourier_n), None)#t_pos_dir)
	def get_P_field_pos(self,n):
		"""
		Gets the sum of all fields projected from the electrodes in the positive direction for fourier compoent n

		Parameters
		----------
		n : int
			fourier compoent

		Returns
		-------
		float
			multiplication factor of projected fields with respect to the case with no reflections
			"""

		if self.P_field_pos[n]==None:
			if not self.P_field_pos[n-1]==1: # if self.P_field_pos[n-1]==1, we conclude that it will be 1 for all further components
				# sum of fields originating from the field initially projected in the positive direction
				if len(self.interfaces)>self.interface_of_electrodes+1:
					rpo=self.interfaces[self.interface_of_electrodes+1].get_r_pos_dir_eff(n)
					rne=self.interfaces[self.interface_of_electrodes].get_r_neg_dir_eff(n)
					t=self.layers[self.interface_of_electrodes+1].t_eff
					self.P_field_pos[n]=1/(1-rpo*rne*np.exp(-np.pi*2*t*(2*n+1)))
				else:
					self.P_field_pos[n]=1
				# sum of fields originating from the field initially projected in the negative direction
				if self.interface_of_electrodes>0:
					rpo=self.interfaces[self.interface_of_electrodes].get_r_pos_dir_eff(n)
					tpo=self.interfaces[self.interface_of_electrodes].get_t_pos_dir_eff(n)
					rne=self.interfaces[self.interface_of_electrodes-1].get_r_neg_dir_eff(n)
					t=self.layers[self.interface_of_electrodes].t_eff
					self.P_field_pos[n]+=rne*tpo*np.exp(-np.pi*2*t*(2*n+1))  /  (1-rpo*rne*np.exp(-np.pi*2*t*(2*n+1)))
				if abs(self.P_field_pos[n]-1)<self.accuracy_limit: # mark that we are done calculating components, by setting self.P_field_pos[n]=1
					self.P_field_pos[n]=1
			else:
				self.P_field_pos[n]=1
		return self.P_field_pos[n]
	def get_P_field_neg(self,n):
		"""
		Gets the sum of all fields projected from the electrodes in the negative direction for fourier compoent n

		Parameters
		----------
		n : int
			fourier compoent

		Returns
		-------
		float
			multiplication factor of projected fields with respect to the case with no reflections
			"""
		if self.P_field_neg[n]==None:
			if not self.P_field_neg[n-1]==1: # if self.P_field_neg[n-1]==1, we conclude that it will be 1 for all further components
				# sum of fields originating from the field initially projected in the positive direction
				if self.interface_of_electrodes>0:
					rpo=self.interfaces[self.interface_of_electrodes].get_r_pos_dir_eff(n)
					rne=self.interfaces[self.interface_of_electrodes-1].get_r_neg_dir_eff(n)
					t=self.layers[self.interface_of_electrodes].t_eff
					self.P_field_neg[n]=1/(1-rpo*rne*np.exp(-np.pi*2*t*(2*n+1)))
				else:
					self.P_field_neg[n]=1
				# sum of fields originating from the field initially projected in the positive direction
				if len(self.interfaces)>self.interface_of_electrodes+1:
					rpo=self.interfaces[self.interface_of_electrodes+1].get_r_pos_dir_eff(n)
					rne=self.interfaces[self.interface_of_electrodes].get_r_neg_dir_eff(n)
					tne=self.interfaces[self.interface_of_electrodes].get_t_neg_dir_eff(n)
					t=self.layers[self.interface_of_electrodes+1].t_eff
					self.P_field_neg[n]+=rpo*tne*np.exp(-np.pi*2*t*(2*n+1))  /  (1-rpo*rne*np.exp(-np.pi*2*t*(2*n+1)))
				if abs(self.P_field_neg[n]-1)<self.accuracy_limit: # mark that we are done calculating components, by setting self.P_field_pos[n]=1
					self.P_field_neg[n]=1
			else:
				self.P_field_neg[n]=1
		return self.P_field_neg[n]
	def get_C(self):
		"""
		Returns
		-------
		float
			capacitace
			"""
		eps_above=self.layers[self.interface_of_electrodes].eps_y*self.layers[self.interface_of_electrodes].eps_r
		eps_below=self.layers[self.interface_of_electrodes+1].eps_y*self.layers[self.interface_of_electrodes+1].eps_r
		return (eps_above+eps_below)*self.Kk/self.Kpk*eps0/2
	def get_C_int_Ex(self):
		"""
		Function for calculating the capacitance by integrating :math:`\\varepsilon_xE_x` at :math:`x=0`.\n
		Used for testing, as it should give same output as get_C().\n
		For all practical applications get_C() should be used instead\n

		Returns
		-------
		float
			capacitace
			"""
		if self.layers[0].eps_x>0:
			G,error=scipy.integrate.quad(lambda y: self.get_Ex(0,y), -np.inf, 0)
			C= G*self.layers[0].eps_x*eps0
		else:
			C=0
		for i in range(self.number_of_layers-2):
			if self.layers[i+1].eps_x>0:
				G,error=scipy.integrate.quad(lambda y: self.get_Ex(0,y), self.interfaces[i].y, self.interfaces[i+1].y)
				C+= G*self.layers[i+1].eps_x*eps0
		if self.layers[-1].eps_x>0:
			G,error=scipy.integrate.quad(lambda y: self.get_Ex(0,y), self.interfaces[-1].y, np.inf)
			C+= G*self.layers[-1].eps_x*eps0
		return C

	def getA(self,n):
		"""
		calculates  :math:`A_{2n+1}`

		Parameters
		----------
		n : fourier compoent

		Returns
		-------
		float
			A_(2n+1)
			"""
		if not n in self.adict:
			Pn=scipy.special.legendre(n)(2*self.k**2-1)
			self.adict[n]=np.pi/self.Kpk*Pn #=A_{2n+1}
		return self.adict[n]


	def get_V_Ex_Ey(self,x,y,get_V=1,get_Ex=1,get_Ey=1): # accepts 'x' as a list, but 'y' must be single value
		"""
		Function for calculating the the potential and electric fields at coordinates (x,y)
		The potential and fields are calculated as the sum of the field going away from the electrodes,
		and the field going towards the electrodes.
		The amplitude of the field going away from the electrode is:
		sum_{n=0->inf} (all t_{eff,n} from electrodes to layer)*exp(-pi*exp_arg_fac*(2n+1))
		where exp(-pi*exp_arg_fac*(2n+1))is calculated as cur_exp=exp(-pi*exp_arg_fac)*exp(-2*pi*exp_arg_fac)^n

		Similarly, the amplitude of the field going away from the electrode is:
		sum_{n=0->inf} (all t_{eff,n} from electrodes to layer)*(r_{eff,n} at the next interface)*exp(-pi*reverse_exp_arg_fac*(2n+1))
		where exp(-pi*reverse_exp_arg_fac*(2n+1))is calculated as cur_exp=exp(-pi*reverse_exp_arg_fac)exp(-2*pi*reverse_exp_arg_fac)^n

		Parameters
		----------
		x : float or list of floats
			x-coordiate(s)
		y : float
			y-coordiate
		get_V : bool, optional
			V is only calculated if this flag is set to True, default: True
		get_Ex : bool, optional
			Ex is only calculated if this flag is set to True, default: True
		get_Ey : bool, optional
			Ey is only calculated if this flag is set to True, default: True

		Returns
		-------
		list of floats for V, Ex, Ey
			"""
		x=np.array(x)
		x=np.atleast_1d(x)
		if self.hybrid:
			if y==self.interfaces[self.interface_of_electrodes].y and get_V==1 and get_Ex==0 and get_Ey==0:
				V=self.get_V_electrodes(x)
				return V, 0,0
		V=np.zeros(x.size)
		Ex=np.zeros(x.size)
		Ey=np.zeros(x.size)
		# get current layer we are in
		layer=0
		while layer<len(self.interfaces) and self.interfaces[layer].y <= y:
			layer+=1
		distPrior=0
		reflected_field=0
		if self.interfaces[self.interface_of_electrodes].y <= y: # we are above electrodes
			direction=1
			getPfield=self.get_P_field_pos
			for i in range(self.interface_of_electrodes+1,layer):
				distPrior+=self.layers[i].t_eff
			# calculate exp_arg_fac, the argument for the exponential, i.e. exp((2n+1)*exp_arg_fac)
			exp_arg_fac=-np.pi*(  (y-self.interfaces[layer-1].y)*self.layers[layer].eps_r+distPrior  )
			if len(self.interfaces) > layer:
				reverse_exp_arg_fac=-np.pi*(  (self.interfaces[layer].y-y)*self.layers[layer].eps_r+self.layers[layer].t_eff+distPrior  )
				reflected_field=1
		else: # we are below electrodes
			direction=-1
			getPfield=self.get_P_field_neg
			for i in range(self.interface_of_electrodes,layer,direction):
				distPrior+=self.layers[i].t_eff
			exp_arg_fac=-np.pi*(  direction*(y-self.interfaces[layer].y)*self.layers[layer].eps_r+distPrior  )
			if layer>0:
				reverse_exp_arg_fac=-np.pi*(  direction*(self.interfaces[layer-1].y-y)*self.layers[layer].eps_r+self.layers[layer].t_eff+distPrior  )
				reflected_field=1
		exp_mul_fac=np.exp(exp_arg_fac*2) # <- multiplication factor for the exponential attenuation
		cur_exp=np.exp(exp_arg_fac) # <- current exponential attenuation
		if reflected_field==1: #if there is a reflected field, do also for this
			reverse_exp_mul_fac=np.exp(reverse_exp_arg_fac*2)
			reverse_cur_exp=np.exp(reverse_exp_arg_fac)
		for n in range(self.max_fourier_n):
			N=2*n+1
			# get the projected field at the interface with the electrodes
			amplitude=self.getA(n)*getPfield(n)
			if direction==1:
				for i in range(self.interface_of_electrodes+1,layer):
					amplitude*=self.interfaces[i].get_t_pos_dir_eff(n)
			else:
				for i in range(self.interface_of_electrodes-1,layer-1,direction):
					amplitude*=self.interfaces[i].get_t_neg_dir_eff(n)
			forward_amplitude=cur_exp*amplitude
			if get_V or get_Ey:
				sin=np.sin(np.pi*N*x)
				if get_V:
					V+=-1/np.pi/N*sin*forward_amplitude
				if get_Ey:
					Ey+=direction*self.layers[layer].eps_r*sin*forward_amplitude
			if get_Ex:
				cos=np.cos(np.pi*N*x)
				Ex+=cos*forward_amplitude
			cur_exp*=exp_mul_fac
			# calculate the reflected field
			reflections_coeff=0
			if direction == 1:
				if len(self.interfaces) > layer:
					reflections_coeff=self.interfaces[layer].get_r_pos_dir_eff(n)
			else:
				if layer>0:
					reflections_coeff=self.interfaces[layer-1].get_r_neg_dir_eff(n)
			if not reflections_coeff==0:
				reverse_amplitude=reverse_cur_exp*reflections_coeff*amplitude
					#sin=np.sin(np.pi*N*x) # <- have already calculated this
				if get_V:
					V+=-1/np.pi/N*sin*reverse_amplitude
				if get_Ey:
					Ey-=direction*self.layers[layer].eps_r*sin*reverse_amplitude
				if get_Ex:
					#cos=np.cos(np.pi*N*x) # <- have already calculated this
					Ex+=cos*reverse_amplitude
				reverse_cur_exp*=reverse_exp_mul_fac
			# if the exponential term is less than accuracy_limit, we are done
			if cur_exp/N<self.accuracy_limit:
				break
		return V, Ex, Ey


	def get_V_inf(self,x):
		"""
		Calculates the potential at the interface with the electrodes
		not considering reflections. ( :math:`V(x,y)` at the
		interface with the electrodes)

		Parameters
		----------
		x : float
			x-coordiate

		Returns
		-------
		float
		 	potential

		"""
		from infinite_conformal import fF
		# offset the geometry to the geometry used by conformal mapping
		x=x+0.5
		# from symmetry we only need 0<x<0.5
		# the period of the system is 2
		x=x%2 # reduce to period
		if x>1:
			x=x%1 # reduce to half the period
			sign_1=-1 # we are on the 'negative' side of half period
		else:
			sign_1=1
		if x>0.5:
			x=1-x # fold onto  0<x<0.5
			sign_2=-1 # we are on the negative side of the fold
		else:
			sign_2=1
		if not x in self.Vdict: # see if we have calculated this before
			z=x+1j*0 # posistion on the z-plane
			t=1/self.k*np.sin(np.pi*z) # posistion on the t-plane
			F=fF(np.arcsin(t),self.k,self.accuracy_limit) # heavy lifting for transformation to w-plane
			self.Vdict[x]=(1-(1/self.Kpk)*F.imag)*0.5
		return self.Vdict[x]*sign_1*sign_2

	def get_V_electrodes(self,x): # accepts 'x' as a list
		"""
		Function for calculating the the potential and electric fields at (x,y)
		at the interface with the electrodes using hybrid formualtion

		Parameters
		----------
		x : 1d numpy array of floats
			x-coordiate(s)

		Returns
		-------
		list of floats for V
			"""

		V=np.zeros(x.size)
		# add initial field using conformal mapping
		for i, XX in enumerate(x):
			V[i]+=self.get_V_inf(XX)
		layer=self.interface_of_electrodes+1
		#z=np.sum([o.z for o in self.interfaces[:layer]])
		# we are above electrodes
		getPfield=self.get_P_field_pos
		if len(self.interfaces) > layer:
			reverse_exp_arg_fac=-np.pi*2*self.layers[layer].t_eff
			reverse_exp_mul_fac=np.exp(reverse_exp_arg_fac*2)
			reverse_cur_exp=np.exp(reverse_exp_arg_fac)
		else:
			reverse_cur_exp=0
		for n in range(self.max_fourier_n):
			N=2*n+1
			P=getPfield(n)
			A2np1=self.getA(n)
			sin=np.sin(N*np.pi*x)
			V+=-1/np.pi/N*A2np1*sin*(P-1)
			R=0
			if len(self.interfaces) > layer:
				R=self.interfaces[layer].get_r_pos_dir_eff(n) #pos dir
			if not R==0:
				V+=-1/np.pi/N*A2np1*sin*reverse_cur_exp*P*R
				reverse_cur_exp*=reverse_exp_mul_fac
			if abs(P-1)/N<self.accuracy_limit and abs(P*reverse_cur_exp/N)<self.accuracy_limit :#np.abs(P-1)/N<self.accuracy_limit:
				break
		return V

	def get_V(self,x,y):
		"""
		calls get_V_Ex_Ey() to calculate V
		"""
		V,Ex,Ey=self.get_V_Ex_Ey(x,y,1,0,0)
		return V
	def get_Ex(self,x,y):
		"""
		calls get_V_Ex_Ey() to calculate Ex
		"""
		V,Ex,Ey=self.get_V_Ex_Ey(x,y,0,1,0)
		return Ex
	def get_Ey(self,x,y):
		"""
		calls get_V_Ex_Ey() to calculate Ey
		"""
		V,Ex,Ey=self.get_V_Ex_Ey(x,y,0,0,1)
		return Ey


class multiple_recursive_images:
	"""
	A class generates and houses instances of the class single_recursive_images and uses a linear combination of these
	to approximate physcal values for the capacitance, potential, and electric fields

	Parameters
	----------
	etas : list of float
		cover fraction of the electrode(s) must be a list of the same length as there are interfaces, i.e. 1 shorter than eps_x_of_layers
		each value must be 0=<eta=<1. 0 indicates that no electrode exists at the corresponding interface
	thickness_of_layers : list of floats
		thicknesses of the layers, this list will be 2 shorter than eps_x_of_layers and eps_y_of_layers, as the outermost layers have no defined thickness, but are infinite
	eps_x_of_layers : list of floats
		in-plane dielectric constant of the layers
	eps_y_of_layers : list of floats
		out-of-plane dielectric constant of the layers
	LAcomp: int or list of ints, optional
		the number of components used for each set of electrodes. If list, the list must be the same length as etas, and have a number larger than one for each nonzero value in etas.
		If not list the same value vill be used for all electrodes. Defaults to 8
	max_fourier_n : int, optional
		maximum number of fourier components. Note that the floating point limit is reached around 190 components, defaults to 120
	voltages : list of floats or None, optional
		If a list of floats is provided, the voltage at each set of electrodes will be set to +/- the corresponding value in the list. the list must have the same length as etas.
		If None, the potential at each set of electrodes with a cover fraction less than 1 is set to +/-0.5, for a total difference of 1V. The potential at electrodes with a cover fraction of 1 is set to 0, as this is assumed to be grounded.
		Defaults to None. NB! the calcualted capacitance assumed a potential of +/-0.5 at the electrodes, and if a different value is used, calculated capacitance will reflect the charge at the electrodes rather than the capacitance!
	accuracy_limit : float, optional
		reflections with amplitude less than accuracy_limit are ignored, defaults to 10**-15

	"""
	def __init__(self,etas,t,eps_x_of_layers,eps_y_of_layers,LAcomp=8,max_fourier_n=120, voltages = None,accuracy_limit=10**-10, hybrid=True):
		# eta = vector of floats, t = vector of floats, eps_x_of_layers = vector of floats, eps_y_of_layers = vector of floats, LAcomp = int, max_fourier_n = int
		# eps_x_of_layers and eps_y_of_layers refer to materials, the length must therefore be at least 2
		# eta refers to interfaces, the length must therefore be at least 1
		# t refers to the thicknes of layers of finite thickness. This vector may have 0 elements.
		# LAcomp must be an int, or a vector of same length as eta
		# voltages must be 'None' or a vector of same length as eta, and is uset to set voltage ratios if multiple sets of electrodes are used
		if not len(eps_x_of_layers)==len(eps_y_of_layers):
			raise NameError('InputSizeError: "eps_x_of_layers" not same size as "eps_y_of_layers"')
		if not len(t)==len(eps_x_of_layers)-2:
			raise NameError('InputSizeError: "t" not 2 shorter than as "eps_y_of_layers"')
		if not len(t)==len(etas)-1:
			raise NameError('InputSizeError: "t" not 1 shortert than "etas"')
		self.accuracy_limit=accuracy_limit
		self.etas=np.array(etas)
		#self.a=1-etas[0]
		#self.b=etas[0]
		self.t=np.array(t)
		self.eps_x_of_layers=np.array(eps_x_of_layers)
		self.eps_y_of_layers=np.array(eps_y_of_layers)
		self.LAcomp=LAcomp
		self.voltages=voltages
		self.max_fourier_n=max_fourier_n
		self.hybrid=hybrid

		# The remainder of the __init__ function is shared with infinite_conformal.multiple_recursive_images.__init()
		#     but is hard-copied in each to increase readability of the code
		if self.voltages == None: # create voltages vector, all electrode sets should have V=0.5, unless they are continious, then they should have 0. If no electrodes are present at the interface, then the voltage is also set to 0
			self.voltages = []
			for eta in self.etas:
				if eta == 0: # no electrodes
					self.voltages.append(0)
				elif eta == 1: # full coverage of electrodes, must have 0 potential
					self.voltages.append(0)
				else:
					self.voltages.append(0.5)
		# declare the step size for for eta in each set of electrodes
		self.electrode_steps=[]
		if isinstance(LAcomp, (int,)):
			for eta in self.etas:
				self.electrode_steps.append(eta/self.LAcomp)
		else:
			for eta, LA in zip(self.etas, self.LAcomp):
				if LA>0:
					self.electrode_steps.append(eta/LA)
				else:
					self.electrode_steps.append(0)
		# calculate the relevant locations and etas for the linear algebra aproach
		self.xpoints=[] #<- list of list of points of interest
		self.xpointVs=[] #<- list of voltages for the linalg expression
		self.CaseEtas=[] #<- list of lists of eta of individual cases
		for eta, step, V in zip(self.etas, self.electrode_steps,self.voltages):
			self.xpoints.append([])
			self.CaseEtas.append([])
			l=step*0.5
			while l<eta and l<1-step: # the second condition prevents the last from beeing added if the electrodes are continious.
				self.xpoints[-1].append(l/2.0-0.5)
				self.CaseEtas[-1].append(l+step*0.5)
				self.xpointVs.append(V)
				l+=step
		self.C=-1
	def get_C(self):
		"""
		Function for calculate the total capacitance using linear algebra.
		To do this it needs to first build the different single_recursive_images objects
		Then use them to calculate the potential at self.xpoints, and use linear algebra to find the weights.
		Finally, the capacitance is calculated by sum(weights*single_recursive_images.get_C)

		Returns
		-------
		float
			capacitace
			"""
		if self.C==-1:
			self.Vs=[]
			self.single=[]
			self.Cs=[]
			# electrodes may exist on multiple interfaces, so iterate
			for interface, _ in enumerate(self.CaseEtas):
				# make the cases with electrodes on this interface
				for  caseEta in self.CaseEtas[interface]:
					if len(self.single)==0:
						self.single.append(single_recursive_images(caseEta,interface,
							self.t,self.eps_x_of_layers,self.eps_y_of_layers,self.max_fourier_n,self.accuracy_limit,
							hybrid=self.hybrid))
					else:
						# utilize the fact that a lot of properties do not need to be recalculated by inheriting from the first case
						self.single.append(single_recursive_images(caseEta,interface,
							self.t,self.eps_x_of_layers,self.eps_y_of_layers,self.max_fourier_n,self.accuracy_limit,
							inherit=self.single[0],hybrid=self.hybrid))
					self.Vs.append(np.array([]))
					# calculate the potential for the last case at all designated locations
					# iterate over designated locations: first iterate over layers
					y=0
					for interfaceOfXpoint, _ in enumerate(self.CaseEtas):
						if len(self.xpoints[interfaceOfXpoint])>0:
							# utilize the fact that get_V suppoprts a list of x coordinates as unput to
							#    calculate all points on this interface by one function call
							self.Vs[-1]=np.concatenate((self.Vs[-1], self.single[-1].get_V(self.xpoints[interfaceOfXpoint],y)))
						if interfaceOfXpoint<len(self.t):
							# need y to keep track of location as we iterate over layers
							y+=self.t[interfaceOfXpoint]
					self.Cs.append(self.single[-1].get_C())
			# solve the linear algebra expression for X: 	AX=B
			A=np.array(self.Vs).transpose()
			B=np.array(self.xpointVs)
			self.weights=numpy.linalg.solve(A,B)
			# calculate the capacitance
			self.C=np.dot(np.array(self.Cs),np.array(self.weights))
		return self.C

	def get_C_int_Ex(self):
		"""
		Function for calculating the capacitance by integrating :math:`\\varepsilon_xE_x` at :math:`x=0`.\n
		Used for testing, as it should give same output as get_C().\n
		For all practical applications get_C() should be used instead\n

		Returns
		-------
		float
			capacitace
			"""
		if self.C==-1:
			self.get_C() # to make sure weights are declared
		Cs=[]
		for case in self.single:
			Cs.append(case.get_C_int_Ex())
		C=np.dot(np.array(Cs),np.array(self.weights))
		return C

	def get_V_Ex_Ey(self,x,y,get_V=1,get_Ex=1,get_Ey=1):
		"""
		Function for calculating the the potential and electric fields at coordinates (x,y)

		Parameters
		----------
		x : float or list of floats
			x-coordiate(s)
		y : float
			y-coordiate
		get_V : bool, optional
			V is only calculated if this flag is set to True, default: True
		get_Ex : bool, optional
			Ex is only calculated if this flag is set to True, default: True
		get_Ey : bool, optional
			Ey is only calculated if this flag is set to True, default: True

		Returns
		-------
		list of float for V, Ex, Ey
			"""
		if self.C==-1:
			self.get_C()
		V=0
		Ex=0
		Ey=0
		for case, weight in zip(self.single, self.weights):
			Vi,Exi,Eyi=case.get_V_Ex_Ey(x,y,get_V,get_Ex,get_Ey)
			V+=Vi*weight
			#print(Vi, weight)
			Ex+=Exi*weight
			Ey+=Eyi*weight
		return(V,Ex,Ey)
	def get_V(self,x,y):
		"""
		calls get_V_Ex_Ey() to calculate V
		"""
		V,Ex,Ey=self.get_V_Ex_Ey(x,y,1,0,0)
		return V
	def get_Ex(self,x,y):
		"""
		calls get_V_Ex_Ey() to calculate Ex
		"""
		V,Ex,Ey=self.get_V_Ex_Ey(x,y,0,1,0)
		return Ex
	def get_Ey(self,x,y):
		"""
		calls get_V_Ex_Ey() to calculate Ey
		"""
		V,Ex,Ey=self.get_V_Ex_Ey(x,y,0,0,1)
		return Ey
