# -1,  3,  5                1 20
# -1,  2, -1                2 0
#  0, -1,  2, -1            3 0
#  0,  0, -1,  2, -1        4 0
#             -1,  2, -1    5 0
#              7,  3,  2    6 55

$d=[-1.,-1.,-1.,-1., 3.]
$a=[-1., 2., 2., 2., 2., 2.]
$h=[ 3.,-1.,-1.,-1.,-1.]
$l=7.
$p=5.
$f=[20.,0.,0.,0.,0.,55.]
$x=res_3diag2(d,a,h,l,p,f)=[1,2,3,4,5,6]
/*
$L=fixx+fiyy=0
*/

class class_potencial {
	function class_potencial($nx,$ny,$lx,$ly):
		self.nx=nx #pocet bodu site ve smeru x
		self.ny=ny
		self.dx=lx/(nx-1)
		self.dy=lx/(ny-1)
		self.dx2=self.dx**2
		self.dy2=self.dy**2
		self.lx=lx
		self.ly=ly
		self.fi=[]
		for iy in range(ny):
			self.fi.append([0.]*nx)
			
	def tisk(self,tx='=========',fmt='{:12.8f}'):
		print(tx)
		for iv in range(self.ny-1,-1,-1):
			ss=''
			for iu in range(self.nx):
				ss+=fmt.format(self.fi[iv][iu])
			print(ss)
	def dfix(self,iu,iv):
		if iu==0: x=-10.
		else:  x=10.
		return(x)
		
	def dfiy(self,iu,iv):
		u=iu*self.dx
		if iv==self.ny-1:
			y=10.
		else:
			y=3*math.exp(-(20.*u-10)**2)
		return(y)
		
	def L(self,ix,iy):
		ff=0.
		if ix==0: 
			ff+=2./3./self.dx2*(self.fi[iy][2]-self.fi[iy][1]-self.dx*self.dfix(ix,iy))
		elif ix==self.nx-1:
			ff+=2./3./self.dx2*(self.fi[iy][ix-2]-self.fi[iy][ix-1]+self.dx*self.dfix(ix,iy))
		else:
			ff+=(self.fi[iy][ix+1]-2.*self.fi[iy][ix]+self.fi[iy][ix-1])/self.dx2
			
		if iy==0:
			ff+=2./3./self.dy2*(self.fi[2][ix]-self.fi[1][ix]-self.dy*self.dfiy(ix,iy))
		elif iy==self.ny-1:
			ff+=2./3./self.dy2*(self.fi[iy-2][ix]-self.fi[iy-1][ix]+self.dy*self.dfiy(ix,iy))
		else:
			ff+=(self.fi[iy+1][ix]-2.*self.fi[iy][ix]+self.fi[iy-1][ix])/self.dy2
		return(ff)

	def af_xiy(self,iy,alfa,omega):
		a=[2.+alfa*self.dx2]*self.nx
		d=[-1.]*(self.nx-1)
		h=[-1.]*(self.nx-1)
		f=[0.]*(self.nx)
		a[0]=3.*self.dx2*alfa
		h[0]=-2.
		p=2.
		
		f[0]=3.*self.dx2*omega*alfa*self.L(0,iy)+2.*self.dx*self.dfix(0,iy)
		for i in range(1,self.nx-1):
			f[i]=omega*alfa*self.dx2*self.L(i,iy)
		f[self.nx-1]=3.*omega*alfa*self.dx2*self.L(self.nx-1,iy)-2.*self.dx*self.dfix(self.nx-1,iy)
		a[self.nx-1]=3.*self.dx2*alfa
		d[self.nx-2]=-2.
		l=2.
		#tiskmat(d,a,h,l,p,f)
		z=res_3diag2(d,a,h,l,p,f)
		#print(z)
		return(z)
		
	def af_x(self,alfa,omega):
		zz=[]
		for iy in range(self.ny):
			#print('******* iy=',iy)
			z=self.af_xiy(iy,alfa,omega)
			zz.append(z)
		return(zz)
		
	def af_yix(self,ix,alfa,zz):
		a=[2.+alfa*self.dy2]*self.ny
		d=[-1.]*(self.ny-1)
		h=[-1.]*(self.ny-1)
		f=[0.]*(self.ny)
		a[0]=3.*self.dy2*alfa
		h[0]=-2.
		p=2.
		f[0]=3.*self.dy2*zz[0][ix] #??+2.*self.dx*self.dfix(0,iy)
		for i in range(1,self.ny-1):
			f[i]=self.dy2*zz[i][ix]
		i=self.ny-1
		a[self.ny-1]=3.*self.dy2*alfa
		d[self.ny-2]=-2.
		l=2.
		f[i]=3.*self.dy2*zz[i][ix] #-2.*self.dx*self.dfix(self.nx-1,iy)
		#tiskmat(d,a,h,l,p,f)
		z=res_3diag2(d,a,h,l,p,f)
		#print(z)
		return(z)
		
	def af_y(self,alfa,zz):
		cc=[]
		for ix in range(self.nx):
			#print('******* ix=',ix)
			c=self.af_yix(ix,alfa,zz)
			cc.append(c)
		return(cc)
	
	def af_step(self,alfa,omega):
		zz=self.af_x(alfa,omega)
		#print(zz)
		cc=self.af_y(alfa,zz)
		#print(cc)
		for ix in range(self.nx):
			for iy in range(self.ny):
				self.fi[iy][ix]+=cc[ix][iy]
		#self.tisk('krok alfa='+str(alfa)+' omega='+str(omega))
		
	def af(self,omega,nalfa):
		dd=math.sqrt(self.nx**2+self.ny**2)
		alh=dd/20.
		ald=1./dd
		print(ald,alh)
		lam=math.pow(alh/ald,1./nalfa)
		alf=[]
		for i in range(nalfa):
			alf.append(ald*math.pow(lam,i))
		print('alf=',alf)
		st=1.e200
		err=1000.
		eps=0.01
		while err>eps:
			for al in alf:
				self.af_step(al,omega)
				#self.tisk(str(al),'{:10.4f}')
				self.odecti(self.fi[self.ny//2][self.nx//2])
				self.tisk('{:10.4f}'.format(al),'{:10.4f}')
				sn=self.norma()
				
				err=abs(sn-st)
				print('{:10.5f} {:10.5f} {:10.5f}'.format(al,sn,self.integral()))
				st=sn
		
	def norma(self):
		for iv in range(self.ny-1,-1,-1):
			ss=0.
			for iu in range(self.nx):
				ss+=self.fi[iv][iu]**2
		return(math.sqrt(ss))
		
	def integral(self):
		ss=0.
		for iv in range(self.ny):
			for iu in range(self.nx):
				ok1=1
				ok2=1
				if iu==0 or iu==self.nx-1: ok1=0
				if iv==0 or iv==self.ny-1: ok2=0
				ok=ok1+ok2
				if ok==2: ss+=self.fi[iv][iu]
				elif ok==1: ss+=self.fi[iv][iu]/2.
				else: ss+=self.fi[iv][iu]/4.
		ss*=self.dx*self.dy
		return(ss)
	
	def odecti(self,delta):
		for iv in range(self.ny-1,-1,-1):
			for iu in range(self.nx):
				self.fi[iv][iu]-=delta
				
	def dej_zebra(self):
		zb=[]
		for j in range(self.ny):
			zz=[]
			for i in range(self.nx):
				zz.append([i*self.dx,j*self.dy,self.fi[j][i]])
			zb.append(zz)
		return(zb)
  }//end of class_potencial
def tiskmat(d,a,h,l,p,f,fmt='{:12.8f}'):
	n=len(a)
	fl=int(fmt[2:fmt.index('.')])
	ss=(fmt+fmt+fmt).format(a[0],h[0],p)
	for j in range(n-3): ss+=' '*fl
	ss+='|'+fmt.format(f[0])
	print(ss)
	for i in range(1,n-1):
		ss=''
		for j in range(i-1): ss+=' '*fl
		ss+=(fmt+fmt+fmt).format(d[i-1],a[i],h[i])
		for j in range(n-i-2): ss+=' '*fl
		ss+='|'+fmt.format(f[i])
		print(ss)
	ss=''
	i=n-1
	for j in range(i-2): ss+=' '*fl
	ss+=(fmt+fmt+fmt).format(l,d[i-1],a[i])
	ss+='|'+fmt.format(f[n-1])
	print(ss)
		
def res_3diag2(d,a,h,l,p,f):
	n=len(a)-1
	h[1]-=d[0]/a[0]*p
	for i in range(n-1):
		w=d[i]/a[i]
		#print('.... i , w',i, w, 'd[',i,']',d[i],'a[',i,']',a[i],h[i])
		a[i+1]-=h[i]*w
		f[i+1]-=f[i]*w
	w=l/a[n-2]
	#print(w,a[n-2],d[n-1],f[n],'l=',l)
	d[n-1]-=w*h[n-2]
	f[n]-=w*f[n-2]
	w=d[n-1]/a[n-1]
	#print(w,a[n-1],d[n-1],f[n],f[n-1],h[n-1])
	a[n]-=h[n-1]*w
	f[n]-=w*f[n-1]
	x=[0.]*(n+1)
	x[n]=f[n]/a[n]
	#print('xn=',x[n],'f=',f[n],'a',a[n])
	for i in range(n-1,0,-1):
		x[i]=(f[i]-h[i]*x[i+1])/a[i]
		#print('** ',i,x[i],'f=',f[i], x[i+1],'h=',h[i],'a=',a[i])
	x[0]=(f[0]-h[0]*x[1]-p*x[2])/a[0]
	return(x)
	

	
#tiskmat(d,a,h,l,p,f)
nx=50
ny=50
fi=class_potencial(nx,ny,1.,1.)
#fi.tisk('1')
omega=1.9
#zz=fi.af_x(alfa,omega)
#print(zz)
#cc=fi.af_y(alfa,zz)
#print(cc)
#fi.af_step(alfa,omega)
#fi.af(1.,10)

grx=[]
grn=[]
gri=[]
grh=[]
al=[200.,100.,50.,25.,12.5,6.25,3.125]
snimek=[]
for k in range(10000):
	kk=k%len(al)
	alfa=al[kk]/min(fi.dx,fi.dy)**2
	fi.af_step(alfa,omega)
	fi.odecti(fi.fi[fi.ny//2][fi.nx//2])
	intg=fi.integral()
	nrm=fi.norma()
	hh=fi.fi[fi.ny-1][0]
	grx.append(k)
	gri.append(intg/10.)
	grn.append(nrm)
	grh.append(hh)
	if k%50: snimek.append(fi)
	
	if k%50==0: print('{:4d}{:10.4f}{:10.4f}{:10.4f}'.format(k,intg,nrm,hh))
	#if k%200==0: 
		#print('*'*50)
		#fi.tisk('{:10.2f}'.format(fi.integral()),'{:8.2f}')
snimek.append(fi)

import matplotlib.pyplot as plt
plt.subplot(3,2,1)
plt.plot(grx,gri,'b', label='integral',lw=1)

plt.subplot(3,2,3)
plt.plot(grx,grn,'g', label='norma',lw=1)
plt.plot(grx,grh,'r', label='hh',lw=1)
plt.legend()

plt.subplot(3,2,2)
v=vr.vrstevnice(snimek[-1].dej_zebra())
v.kresli_vrstevnice(5,'black',1)

xx=[]
uu=[]
yy=[]
vv=[]
for iy in range(fi.ny):
	wx=[]
	wy=[]
	wu=[]
	wv=[]
	for ix in range(fi.nx):
		if ix==0: dx=(-3.*fi.fi[iy][0]+4.*fi.fi[iy][1]-fi.fi[iy][2])/2./fi.dx
		elif ix==fi.nx-1: dd=(fi.fi[iy][fi.nx-3]-4.*fi.fi[iy][fi.nx-2]+3.*fi.fi[iy][fi.nx-1])
		else: dx=(fi.fi[iy][ix+1]-fi.fi[iy][ix-1])/2./fi.dx
		
		if iy==0: dy=(-3.*fi.fi[0][ix]+4.*fi.fi[1][ix]-fi.fi[2][ix])/2./fi.dy
		elif iy==fi.ny-1: dy=(fi.fi[fi.ny-3][ix]-4.*fi.fi[fi.ny-2][ix]+ 3.*fi.fi[fi.ny-1][ix])/2./fi.dy
		else: dy=(fi.fi[iy+1][ix]-fi.fi[iy-1][ix])/2./fi.dy
		u=ix*fi.dx
		v=iy*fi.dy
		wx.append([u,v,dx])
		wy.append([u,v,dy])
		wu.append([dx,dy,u])
		wv.append([dx,dy,v])
	xx.append(wx)
	yy.append(wy)
	uu.append(wu)
	vv.append(wv)
	
plt.subplot(3,2,4)
xv=vr.vrstevnice(xx)
xv.kresli_vrstevnice(10,'g',1)
yv=vr.vrstevnice(yy)
yv.kresli_vrstevnice(10,'b',1)

plt.subplot(3,2,6)
xv=vr.vrstevnice(uu)
xv.kresli_vrstevnice(10,'g',1)
yv=vr.vrstevnice(vv)
yv.kresli_vrstevnice(10,'b',1)
yv.kresli_vrstevnice([0.001,9.999],'b',1)

plt.show()