#############################################################################################################
#																											#
#	Title			:	Simulation of wave function using standard volume approach for circular quantum wir-#
#						e geometry mesh generated by Gmsh and visualized									#
#	Author			:	Alfred Ajay Aureate R																#
#	Roll No			:	EE10B052																			#
#	Project guide	:	Prof. Anil Prabhakar, Dept.of Electrical Engineering, IITM							#
#	Code location	:	/DDP_codes/sim_mesh/circular_qw/cqw.py												#
#	Figure ref.		:	NA																					#
#	Date			:	19th May, 2015																		#
#																											#
#############################################################################################################

from fipy import *
import numpy as np
cellSize = 0.05
radius = 1.

#############################################################################################################
# generates circular mesh																					#
#############################################################################################################
mesh = Gmsh2D('''
			cellSize = %(cellSize)g;
			radius = %(radius)g;
			Point(1) = {0, 0, 0, cellSize};
			Point(2) = {-radius, 0, 0, cellSize};
			Point(3) = {0, radius, 0, cellSize};
			Point(4) = {radius, 0, 0, cellSize};
			Point(5) = {0, -radius, 0, cellSize};
			Circle(6) = {2, 1, 3};
			Circle(7) = {3, 1, 4};
			Circle(8) = {4, 1, 5};
			Circle(9) = {5, 1, 2};
			Line Loop(10) = {6, 7, 8, 9};
			Plane Surface(11) = {10};
			'''%locals())
phi = CellVariable(name = "solution variable", mesh = mesh, value = 0.)
D = 1.
eq = ImplicitSourceTerm(coeff=0.21742) == DiffusionTerm(coeff=D)	# don't know why 0.21742 gives normalized solution, but it does!

#############################################################################################################
# setting the boundary conditions (shifted down by 1 later)													#
#############################################################################################################
valueExterior = 1
phi.constrain(valueExterior, mesh.exteriorFaces)

phi1 = 1-phi														# shifted down by 1 for correctness of the expected boundary conditions
phi1.setName("$\psi$")
phis = abs(phi1*phi1).sum()*(np.pi)	
phi2 = phi1/np.sqrt(phis)											# normalizing (temporarily)
phi2.setName("$\psi$")

#############################################################################################################
# Visualizing the ground state of the circular quantum wire for about 5 seconds								#
#############################################################################################################
if __name__ == '__main__':
	viewer = Viewer(vars=phi2, datamin=0., datamax=0.02)
	viewer.plot()
timeStepDuration = 10 * 0.9 * cellSize**2 / (2 * D)
steps = 10
for step in range(steps):
	eq.solve(var=phi, dt=timeStepDuration)
	phi1 = 1-phi
	phis = abs(phi1*phi1).sum()*(np.pi)
	phi2 = phi1/np.sqrt(phis)											# normalizing (temporarily)
	if __name__ == '__main__':
		viewer.plot()

