% --- Executes on button press in capacitance.
function capacitance_Callback(hObject, eventdata, handles)
global electrode mx my CAP analyse_press cap_press LCAP pot_press
tic 
colorbar off
cla
load analys_var ID LOC
pot_press=0;
overlay()
HH= findobj(gcf,'Tag','cap');
set(HH,'String',' ')

[disc,page]=size(electrode);
for i=1:page
    electrode(i).number=i;
end

flag1=0;
flag2=0;
flag3=0;

[xinput,yinput]=ginput(1);
x=round(xinput);
y=round(yinput);
for i=1:page
    pep=electrode(i).elem;
    [mpep,npep]=size(pep);
    for j=1:mpep
        if (pep(j,1)==x)&&(pep(j,2)==y)
            flag1=1;
            a=electrode(i).number;
        end
    end
end

if flag1==0
    msgbox('Error: No electrode at selected point')
else
    [xinput,yinput]=ginput(1);
    x=round(xinput);
    y=round(yinput);
    for i=1:page
        pep=electrode(i).elem;
        [mpep,npep]=size(pep);
        for j=1:mpep
            if (pep(j,1)==x)&&(pep(j,2)==y)
                pterr=0;
                flag2=1;
                b=electrode(i).number;
            end
        end
    end
end

if (flag1==1)&&(flag2==0)
    msgbox('Error: No electrode at selected point')
end

delx=1;    
dely=1;    
beta=delx/dely;               
betasq=beta^2;              
gama=-2*(1+betasq);         
[no page]=size(electrode);

CAP=zeros(page,page);
location=max(max(LOC));

if (analyse_press==1)&&(flag1==1)&&(flag2==1)&&(cap_press==0)
    flag3=1;
    msgbox('Press "OK", then wait till this msgbox disappears. Please wait till all the electrodes are numbered and the answer is displayed in the box')
    for laser=1:page
        BC=zeros(mx,my);
        temp=electrode(laser).elem;
        [m,n]=size(temp);
        for i=1:m
            BC(temp(i,2),temp(i,1))=1;
        end
        maxwell(BC,laser,location,gama,betasq)
        euro=num2str(electrode(laser).number);
        x=temp(1,1);
        y=temp(1,2);
        [xc,yc]=findmid(x,y);
        axes(handles.grid)
        text(xc,yc,euro,'Color','r','FontWeight','bold','Fontsize',12)
    end
elseif analyse_press~=1
    msgbox('Press Analyse first');
end

% for i=1:page
%    V(i,1)=electrode(i).voltage; 
% end
% Q=CAP*V

if flag3==1
    LCAP=zeros(page,page);
    for i=1:page
        for j=1:page
            if j==i
                LCAP(i,j)=sum(CAP(i,:));
            else
                LCAP(i,j)=-CAP(i,j);
            end
        end
    end
    eval=num2str(LCAP(a,b));
%     disp('Maxwell/ground capacitance matrix is:')
%     disp(CAP)
%     disp('Lumped capacitance matrix is:')
%     disp(LCAP)
elseif cap_press==1
    %if 'capacitance' was already pressed once, then just number the
    %electrodes and display the capacitance from the already formed matrix
    [disc,page]=size(electrode);
    for i=1:page
        temp=electrode(i).elem;
        x=temp(1,1);
        y=temp(1,2);
        [xc,yc]=findmid(x,y);
        euro=num2str(electrode(i).number);
        axes(handles.grid)
        text(xc,yc,euro,'Color','w','FontWeight','bold','Fontsize',12)
    end
    eval=num2str(LCAP(a,b));
else
    eval=' ';
end

%HH=findobj(gcf,'Tag','cap')
set(HH,'String',eval)

if cap_press==0
    cap_press=cap_press+1;
end
toc


function maxwell(BC,laser,location,gama,betasq)
global CAP mx my electrode shape ID invisible
load analys_var ID LOC
delx=1;    
dely=1;
COF=zeros(location,location);
FV=zeros(location,1);

%ASSEMBLING CO-EFFICIENT MATRIX FOR CALCULATION OF POTENTIAL DISTRIBUTION%
for i=1:mx
    for j=1:my
        ptr=LOC(i,j);
        if ptr==0
            continue
        elseif (i==1)||(j==1)||(i==mx)||(j==my)
            if (i==1)&&(j>=1)&&(j<=my)
                COF(ptr,ptr)=COF(ptr,ptr)-1;
                COF(ptr,LOC(i+1,j))=1;
            end
            if (i==mx)&&(j>=1)&&(j<=my)
                COF(ptr,ptr)=COF(ptr,ptr)-1;
                COF(ptr,LOC(i-1,j))=1;
            end
            if (i>=1)&&(i<=mx)&&(j==1)
                COF(ptr,ptr)=COF(ptr,ptr)-1;
                COF(ptr,LOC(i,j+1))=1;
            end
            if (i>=1)&&(i<=mx)&&(j==my)
                COF(ptr,ptr)=COF(ptr,ptr)-1;
                COF(ptr,LOC(i,j-1))=1;
            end
        else
            COF(ptr,ptr)=gama;
            if ID(i,j+1)==1
                FV(ptr)=FV(ptr)-((BC(i,j+1))*1);
            else
                COF(ptr,LOC(i,j+1))=1;
            end
            
            if ID(i,j-1)==1
                FV(ptr)=FV(ptr)-((BC(i,j-1))*1);
            else
                COF(ptr,LOC(i,j-1))=1;
            end
        
            if ID(i+1,j)==1
                FV(ptr)=FV(ptr)-((BC(i+1,j))*betasq);
            else
                COF(ptr,LOC(i+1,j))=betasq;
            end
                      
            if ID(i-1,j)==1
                FV(ptr)=FV(ptr)-((BC(i-1,j))*betasq);
            else
                COF(ptr,LOC(i-1,j))=betasq;
            end
        end
    end
end
Q=COF\FV;

%ASSEMBLING POTENTIAL MATRIX%
POT=zeros(mx,my);
for i=1:mx
    for j=1:my
        ptr=LOC(i,j);
        if ptr==0
            POT(i,j)=BC(i,j);
        else
            POT(i,j)=Q(ptr);
        end
    end
end

%ASSEMBLING ELECTRIC FIELD MATRIX, WITH DIRECTION
[DX,DY]=gradient(POT);
DX=-DX;
DY=-DY;
EF=sqrt(DX.^2+DY.^2);
for i=1:mx
    for j=1:my
        if BC(i,j)==0
            EF(i,j)=-EF(i,j);
        end
    end
end

%OBTAINING POINTS ONLY ON ELECTRODE BOUNDARIES AND 
%INTEGRATING THE CHARGE AT THESE POINTS
boundary=[0 0 0];
[disc,page]=size(shape);
epsilon=8.854*1e-12;
for i=1:page
    temp=shape(i).elem;
    col1=temp(:,1);
    col2=temp(:,2);
    x=min(col1);
    y=min(col2);
    xmax=max(col1);
    ymax=max(col2);
    w=xmax-x;
    h=ymax-y;
    if shape(i).id==0           %ELLIPSE
        xc=(x+xmax)/2;
        yc=(y+ymax)/2;
        t=1;
        w=w/2;
        h=h/2;
        for s=x:xmax
            q=(sqrt(1-(((s-xc)/w)^2)))*h;
            a=yc+q;
            ell(t,1)=s;
            ell(t,2)=a;
            t=t+1;
        end
        for s=(xmax-1):-1:(x+1)
            q=(sqrt(1-(((s-xc)/w)^2)))*h;
            b=yc-q;
            ell(t,1)=s;
            ell(t,2)=b;
            t=t+1;
        end
        t=t-1;
        for r=1:t
            if r==1
                x=ell(1,1);
                y=ell(1,2);
                dist=sqrt(((ell(t,1)-ell(1,1))^2)*delx+((ell(t,2)-ell(1,2))^2)*dely);
                EFI=EF(round(y),round(x));
                ell(1,3)=EFI*epsilon*dist*2;
            else
                x=ell(r,1);
                y=ell(r,2);
                dist=sqrt(((ell(r,1)-ell(r-1,1))^2)*delx+((ell(r,2)-ell(r-1,2))^2)*dely);
                EFI=EF(round(y),round(x));
                ell(r,3)=EFI*epsilon*dist*2;
            end
        end
        boundary=vertcat(boundary,ell);
        clear ell
    elseif shape(i).id==1      %RECTANGLE
        t=1;
        for s=x:xmax
           rect(t,1)=s;
           rect(t,2)=ymax;
           t=t+1;
        end
        for s=(ymax-1):-1:y
            rect(t,1)=xmax;
            rect(t,2)=s;
            t=t+1;
        end
        for s=(xmax-1):-1:x
            rect(t,1)=s;
            rect(t,2)=y;
            t=t+1;
        end
        for s=(y+1):(ymax-1)
            rect(t,1)=x;
            rect(t,2)=s;
            t=t+1;
        end
        t=t-1;
        for r=1:t
            if r==1
                x=rect(1,1);
                y=rect(1,2);
                dist=sqrt(((rect(t,1)-rect(1,1))^2)*delx+((rect(t,2)-rect(1,2))^2)*dely);
                EFI=EF(round(y),round(x));
                rect(1,3)=EFI*epsilon*dist*2;
            else
                x=rect(r,1);
                y=rect(r,2);
                dist=sqrt(((rect(r,1)-rect(r-1,1))^2)*delx+((rect(r,2)-rect(r-1,2))^2)*dely);
                EFI=EF(round(y),round(x));
                rect(r,3)=EFI*epsilon*dist*2;
            end
        end
        boundary=vertcat(boundary,rect);
        clear rect 
    end
    [disc,pagei]=size(invisible);
    hut=0;
    for j=1:pagei
        if invisible(j).shape==shape(i).seq
            hut=j;
        end
    end
    if hut~=0
        temp=invisible(hut).elem;
        col1=temp(:,1);
        col2=temp(:,2);
        x=min(col1);
        y=min(col2);
        xmax=max(col1);
        ymax=max(col2);
        w=xmax-x;
        h=ymax-y;
        if invisible(hut).id==0           %ELLIPSE
            xc=(x+xmax)/2;
            yc=(y+ymax)/2;
            t=1;
            w=w/2;
            h=h/2;
            for s=x:xmax
                q=(sqrt(1-(((s-xc)/w)^2)))*h;
                a=yc+q;
                ell(t,1)=s;
                ell(t,2)=a;
                t=t+1;
            end
            for s=(xmax-1):-1:(x+1)
                q=(sqrt(1-(((s-xc)/w)^2)))*h;
                b=yc-q;
                ell(t,1)=s;
                ell(t,2)=b;
                t=t+1;
            end
            t=t-1;
            for r=1:t
                if r==1
                    x=ell(1,1);
                    y=ell(1,2);
                    dist=sqrt(((ell(t,1)-ell(1,1))^2)*delx+((ell(t,2)-ell(1,2))^2)*dely);
                    EFI=EF(round(y),round(x));
                    ell(1,3)=EFI*epsilon*dist*2;
                else
                    x=ell(r,1);
                    y=ell(r,2);
                    dist=sqrt(((ell(r,1)-ell(r-1,1))^2)*delx+((ell(r,2)-ell(r-1,2))^2)*dely);
                    EFI=EF(round(y),round(x));
                    ell(r,3)=EFI*epsilon*dist*2;
                end
            end
            boundary=vertcat(boundary,ell);
            clear ell
        elseif invisible(hut).id==1      %RECTANGLE
            t=1;
            for s=x:xmax
               rect(t,1)=s;
               rect(t,2)=ymax;
               t=t+1;
            end
            for s=(ymax-1):-1:y
                rect(t,1)=xmax;
                rect(t,2)=s;
                t=t+1;
            end
            for s=(xmax-1):-1:x
                rect(t,1)=s;
                rect(t,2)=y;
                t=t+1;
            end
            for s=(y+1):(ymax-1)
                rect(t,1)=x;
                rect(t,2)=s;
                t=t+1;
            end
            t=t-1;
            for r=1:t
                if r==1
                    x=rect(1,1);
                    y=rect(1,2);
                    dist=sqrt(((rect(t,1)-rect(1,1))^2)*delx+((rect(t,2)-rect(1,2))^2)*dely);
                    EFI=EF(round(y),round(x));
                    rect(1,3)=EFI*epsilon*dist*2;
                else
                    x=rect(r,1);
                    y=rect(r,2);
                    dist=sqrt(((rect(r,1)-rect(r-1,1))^2)*delx+((rect(r,2)-rect(r-1,2))^2)*dely);
                    EFI=EF(round(y),round(x));
                    rect(r,3)=EFI*epsilon*dist*2;
                end
            end
            boundary=vertcat(boundary,rect);
            clear rect 
        end
    end
end
boundary(1,:)=[];

%STORE BOUNDARY ELEMENTS OF INDIVIDUAL ELECTRODES IN STRUCTURE
[mfel,nfel]=size(boundary);
for i=1:mfel
    boundary(i,2)=round(boundary(i,2));
end
bounds=boundary;
bounds(:,3)=[];

[disc,page]=size(electrode);
for i=1:page
    temp=electrode(i).elem;
    train=intersect(temp,bounds,'rows');
    electrode(i).boundary=train;
    q=0;
    [mtrain,ntrain]=size(train);
    [mboundary,nboundary]=size(boundary);
    for t=1:mtrain
       x=train(t,1);
       y=train(t,2);
       for p=1:mboundary
           if (boundary(p,1)==x)&&(boundary(p,2)==y)
               q=q+boundary(p,3);
           end
       end
    end
    electrode(i).charge=q;
end

%FORM THE iTH COLUMN OF MAXWELL CAPACITANCE MATRIX
for i=1:page
   Qcharge(i,1)=electrode(i).charge;
end
CAP(1:page,laser)=Qcharge(1:page);
clear Qcharge
return

