
% make the grid 

pot_grid=zeros(100, 100); 

% % SET boundary conditions
% pot_grid(100,:)=200; % bottom boundary (target)
% pot_grid(:,1)=0; % left boundary
% pot_grid(:,100)=0; % right boundary
% pot_grid(1,:)=0; % top boundary
% pot_grid(75,:)=0; % plasma sheath boundary
% 
% % potential hump
% % pot_grid(500,500)=20; 
% % pot_grid(500,490)=15; 
% % pot_grid(500,480)=10; 
% % pot_grid(500,470)=5; 
% 
% pot_grid(1,1)=0; % top left corner
% pot_grid(1,100)=0; % top right corner
% pot_grid(100,1)=100; % bot left corner
% pot_grid(100,100)=100; % bot right corner
figure('units','normalized','outerposition',[0 0 1 1]);
avg_kernel=fspecial('average',3); %[0 1/4 0; 1/4 0 1/4; 0 1/4 0];
new_pot_grid=[];
waiter=waitbar(0, 'wait'); 
steps=10000; 
for t=1:steps
    
% RESET boundary conditions
pot_grid(40,1:20)=0; % top left 
pot_grid(40,30:60)=200; % top middle
pot_grid(40,70:100)=0; % top right
% bottom 
pot_grid(60,1:20)=0; % bottom left 
pot_grid(60,30:60)=200; % bottom middle
pot_grid(60,70:100)=0; % bottom right
% pot_grid(100,:)=0; % bottom boundary (target)
% pot_grid(:,1)=0; % left boundary
% pot_grid(:,100)=0; % right boundary
% pot_grid(1,:)=0; % top boundary
% pot_grid(95,:)=0; % plasma sheath boundary

% potential hump
% pot_grid(50,50)=20; 
% pot_grid(50,45)=15; 
% pot_grid(50,40)=10; 
% pot_grid(50,35)=5; 

% pot_grid(90,50)=20; 
% pot_grid(90,45)=15; 
% pot_grid(90,40)=10; 
% pot_grid(90,35)=5; 

% pot_grid(1,1)=0; % top left corner
% pot_grid(1,100)=0; % top right corner
% pot_grid(100,1)=100; % bot left corner
% pot_grid(100,100)=100; % bot right corner

% average_kernel_2=[0 1/6 0; 1/6 0 1/6; 0 1/6 0];
% average_kernel_1=[0 0 0; 0 1/6 0; 0 0 0];
% average_kernel_3=[0 0 0; 0 1/6 0; 0 0 0];
% average_kernel=cat(3, average_kernel_1, average_kernel_2, average_kernel_3);

new_pot_grid=conv2(pot_grid, avg_kernel); 
temp_pot_grid=new_pot_grid(2:101,2:101);
% abs(max(max(temp_pot_grid-pot_grid)))
% if abs(max(max(temp_pot_grid-pot_grid)))<1
% break
% end
pot_grid=new_pot_grid(2:101,2:101); 

waitbar(t/steps, waiter); 
% contour(1:100, 1:100, flipud(pot_grid)); title(sprintf('new pot grid with steps %f', steps)) 
if mod(t,5)==0 && t<300
contour(pot_grid(1:end,1:end)); shading flat; hold on; title(sprintf('%g',t)); 
axis([0 100 0 100 0 200]); 
saveas(gcf,['potential_',sprintf('timestep equals %g',t),'.tiff']); 
end

end
close(waiter)

pot_grid=flipud(pot_grid); 
% interpolate the pot_grid
[xx yy]=meshgrid(1:0.1:100, 1:0.1:100); 
x=1:100; y=1:100; 
pot_grid=interp2(x, y, pot_grid, xx, yy); 
B_grid=ones(size(pot_grid)); 
figure; contour(pot_grid(1:end,1:end)); shading flat; axis equal; hold on; 

% calculate E field 
[EgridX EgridY]=gradient(pot_grid);
EgridX= -40*EgridX; EgridY= -40*EgridY; 
% figure; contour(xx, yy, pot_grid); hold on; axis equal; 
% quiver(EgridX(1:10:end,1:10:end), EgridY(1:10:end,1:10:end))
% quiver(EgridX, EgridY); 


% contour(1:100, 1:100, pot_grid); title(sprintf('new pot grid with steps %f', steps)) 
% hold on; 
% [delx dely]=gradient(pot_grid, 1, 1); 
% quiver(1:100, 1:100,-delx(2:101), -dely(2:101)); 

% imshow(new_pot_grid); shading flat; colormap(jet)
% imcontrast; 

% set(gcf, 'units', 'normalized', 'outerposition', [0 0 1 1]);   

% % figure; imshow(pot_grid); title('original pot grid')
% imshow(new_pot_grid); title(sprintf('new pot grid with steps %f', steps)) 
% colormap(jet)
% impixelinfo; 
% imcontrast; 
% 
% figure; imshow(-1E5*del2(new_pot_grid(10:90, 10:90), 3)); title(sprintf('E field with steps %f', steps))
% imcontrast; 
% colormap(jet)
% impixelinfo; 
% 




