function nma_MAE185_proj_4()
%
% SOlves project #4 for course MAE 185
% UCI, spring 2003.
%
% This solves the second problem in project #4.
%
% See report for full description of problem and
% output.
%
% AUthor: Nasser Abbasi


FALSE = 0;
TRUE  = 1;

TOLERANCE=0.0001;

% since we solve for one quadrant only
% 10 cm x 10 cm.
%
nRow = 11;
nCol = 11;

Z = zeros(nRow,nCol);
W = zeros(nRow,nCol);

% Now add the voltage at the correct location in the quadrant.

Z(nRow-1,1:3) = 100;
Z(nRow-2,1:3) = 100;

%Now solve for the Z plane

while(1)
    Z_old=Z;            
    
    for(i=1:nRow-1)
        for(j=1:nCol)
            if( (i == 10 | i==9)  & (j==1 | j==2 | j==3))
                ;
            else
                if(i==1)
                    top = W(end-1,j);
                    if(j~=nCol)
                        right=Z(i+1,j);
                    else
                        right = W(end-i+1,end-1);
                    end
                    btm=Z(i+1,j);
                    if(j~=1)
                        left=Z(i,j-1);
                    else
                        left=Z(i,j+1);
                    end
                else
                    if(j==1)
                        left = Z(i,j+1);
                        right=Z(i,j+1);
                        top=Z(i-1,j);
                        btm=Z(i+1,j);
                    else
                        if(j==nCol)
                            right = W(end-i+1,end-1);
                            top=Z(i-1,j);
                            left=Z(i,j-1);
                            btm=Z(i+1,j);
                        else
                            top=Z(i-1,j);
                            left=Z(i,j-1);
                            right=Z(i,j+1);
                            btm=Z(i+1,j);
                        end
                    end                                                            
                end
                
                Z(i,j)=top+btm+left+right;
                Z(i,j)=Z(i,j)/4;                
            end                        
        end
    end
    
    C=Z;
    C(nRow-1,1:3) = 0;
    C(nRow-2,1:3) = 0;    
    [maxRow,maxRowIndex]=max(C);
    [maxEle,maxColIndex]=max(maxRow);
    pos=[maxRowIndex(1),maxColIndex(1)];
    x=pos(1);
    y=pos(2);
    eps=abs(   (Z(x,y) - Z_old(x,y)) / Z(x,y) ) * 100;
    if( eps < TOLERANCE )
        break;
    end
end

% Now we solve for the W plane. The edges of the W plane
% are the same as the edges of the Z plane, but put in the
% correct ordering when copying.
%

W(end,:) = Z(1,:);

for(i=1:nRow)
    W(end-(i-1),end)=Z(i,end);
end


%Ok, edges of the W plane are in, now solve for 
%the inside of the W plane

while(1)
    W_old=W;            
    
    for(i=2:nRow-1)
        for(j=1:nCol-1)
            if(j==1)
                left = W(i,j+1);
                right=W(i,j+1);
                top=W(i-1,j);
                btm=W(i+1,j);
            else
                top=W(i-1,j);
                left=W(i,j-1);
                right=W(i,j+1);
                btm=W(i+1,j);
            end
            if(i==9 & j==9)
                i;
            end
            
            W(i,j)=top+btm+left+right;
            W(i,j)=W(i,j)/4;                
        end
    end
    
    [maxRow,maxRowIndex]=max(W(1:end-1,1:end-1));
    [maxEle,maxColIndex]=max(maxRow);
    pos=[maxRowIndex(1),maxColIndex(1)];
    x=pos(1);
    y=pos(2);
    eps=abs(   (W(x,y) - W_old(x,y)) / W(x,y) ) * 100;
    if( eps < TOLERANCE)
        break;
    end
end

%NOw W and Z are solves. To get the full solution
%Now add the other halfs from symmetry

Zfull = zeros( nRow,(nCol+nCol-1) );
Wfull = zeros( nRow,(nCol+nCol-1) );

for(i=1:nRow)
    k=0;
    for(j=nCol:-1:1)
        k=k+1;
        Zfull(i,j)=Z(i,k);
        Wfull(i,j)=W(i,k);
    end
end
Zfull(:,nCol:end)=Z(:,:);
Wfull(:,nCol:end)=W(:,:);


% Now we have full solution completed in Zfull and Wfull.
% The rest is just plotting

figure;
subplot(2,1,1);
mesh(Zfull)
view(-116,64)
title('Z plane');
xlabel('x');
ylabel('y');
zlabel('Voltage');
%setAxisLimits;

subplot(2,1,2);
mesh(Wfull)
view(-109,52)
title('W plane');
xlabel('x');
ylabel('y');
zlabel('Voltage');
%setAxisLimits;


plotit(Zfull,'Z');
plotit(Wfull,'W');

figure;
subplot(2,1,1);
[CS,H]=contour(Zfull);
clabel(CS,H,'FontSize',8); % using H makes the label show inside the lines
title('Contour plot of the Voltage level in the Z plane');
setAxisLimits;
xlabel('x cm');
ylabel('y cm');
colorbar;

subplot(2,1,2);
[CS,H]=contour(Wfull);
clabel(CS,H,'FontSize',8);
title('Contour plot of the Voltage level in the W plane');
setAxisLimits;
xlabel('x cm');
ylabel('y cm');
colorbar;


figure;
subplot(2,1,1);
contourf(Zfull); 
title('Full Contour plot of the Voltage level in the Z plane');
setAxisLimits;
xlabel('x cm');
ylabel('y cm');
colorbar;

subplot(2,1,2);
contourf(Wfull); 
title('Full Contour plot of the Voltage level in the W plane');
setAxisLimits;
xlabel('x cm');
ylabel('y cm');
colorbar;


%%%%%%%%%%%%%%
%
%
%%%%%%%%%%%%%%%
function plotit(M,name)

figure;
[nRowsFull,nColsFull]=size(M);
plot(nRowsFull,nColsFull,'*');
hold on;
plot(1,1,'*');
XLim([1,nColsFull]);
YLim([1,nRowsFull]);

for(i=1:nRowsFull)
    yCoord=nRowsFull-i+1;
    for(j=1:nColsFull)
        xCoord= nColsFull-j+1;
        val=sprintf('%.3f',M(i,j));
        text(xCoord,yCoord,val,'FontSize',7);    
    end
end

title(sprintf('%s plane Voltage values',name));
xlabel('X cm');
ylabel('Y cm');

setAxisLimits;

%%%%%%%%%%%%%%%%%%
%
%
%%%%%%%%%%%%%%%%%%%%5
function setAxisLimits()

newXlabel=[1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19];
l=get(gca,'xtick');
l=linspace(l(1),l(end),length(newXlabel));
set(gca,'xtick',l);
set(gca,'xticklabel',newXlabel);

newYlabel=[0 1 2 3 4 5 6 7 8 9 10];
l=get(gca,'ytick');
l=linspace(l(1),l(end),length(newYlabel));
set(gca,'ytick',l);
set(gca,'yticklabel',newYlabel);