%%%%%%%%%%%%%%%%%%%%%%%%%%
%   vorgumeet1.m
%

clear
a=10; %m=100; 
% im=1;
k1wv=[]; mv=[];
    rajat=1,    % 1-jäik    2-vaba
    kgr=1
% for m=[4:2:20,24:4:60,70:10:100]
for m=24
    n=m
    M=m+2;  N=n+2;
    dx=a/m; dy=dx;  al=dx/dy;
    xv=0:dx:a+.00001;
    yv=0:dy:a+.00001;
    XV=[-dx,xv,a+dx];
    YV=[-dy,yv,a+dy];
    [xm,ym]=meshgrid(xv,yv);
    
    if kgr==1
        figure(1)
        plot([xv(1), xv(end), xv(end), xv(1), xv(1)],...
            [yv(1),  yv(1),  yv(end),  yv(end),yv(1)],'b-','linewidth',3)
        hold on
        plot(xm,ym,'r.','markersize',20)
        plot(xv(2:end-1),ym(1,2:end-1)-dy,'ro')
        plot(xv(2:end-1),ym(end,2:end-1)+dy,'ro')
        plot(xm(2:end-1,1)-dx,yv(2:end-1),'ro')
        plot(xm(2:end-1,end)+dx,yv(2:end-1),'ro')
        plot([xm(2:end-1,1)-dx,xm(2:end-1,end)+dx]',[yv(2:end-1)',yv(2:end-1)']','r--')
        plot([xv(2:end-1)',xv(2:end-1)']',[ym(1,2:end-1)-dy;ym(end,2:end-1)+dy],'r--')
        plot([xv(1),xv(1)]-dx,[yv(2),yv(end-1)],'r--')
        plot([xv(end),xv(end)]+dx,[yv(2),yv(end-1)],'r--')
        plot([xv(2),xv(end-1)],[yv(1),yv(1)]-dy,'r--')
        plot([xv(2),xv(end-1)],[yv(end),yv(end)]+dy,'r--')
        plot([xm(2:end-1,1),xm(2:end-1,end)]',[yv(2:end-1)',yv(2:end-1)']','r-')
        plot([xv(2:end-1)',xv(2:end-1)']',[ym(1,2:end-1);ym(end,2:end-1)],'r-')
        axis([xv(1)-1.1*dx,xv(end)+1.1*dx,yv(1)-1.2*dy,yv(end)+1.3*dy])
        axis equal
        axis ij
        hold off
    end
    % grid
    
    textm=zeros(M+1,N+1);
    i=N/2+1;  j=M/2+1; textm(i,j)=1;
    tii=textm(i,i);
    NN=1;
    for i=M/2 : -1 : 1
        %     i
        jj=M/2+1:-1:M/2-NN+1;
        textm(i,jj)=tii+1:tii+length(jj);
        textm(jj,i)=(tii+1:tii+length(jj))';
        tii=textm(i,i);
        %     i=i-1
        NN=NN+1;
    end
    for j=1:N/2+1
        textm(M+1:-1:M/2+2,j)=textm(1:N/2,j);
    end
    for i=1:M+1
        textm(i,N+1:-1:N/2+2)=textm(i,1:M/2);
    end
    textm(1,[1 2 N N+1])=[0 0 0 0];
    textm(M+1,[1 2 N N+1])=[0 0 0 0];
    textm([1 2 N N+1],1)=[0 0 0 0]';
    textm([1 2 N N+1],N+1)=[0 0 0 0]';
    % textm
    if kgr==1
        for i=1:N+1
            for j=1:M+1
                if textm(i,j)~=0, h=text(XV(i)+.1*dx,YV(j)+.15*dy,int2str(textm(i,j))); end
            end
        end
        axis off
        times12
        print -depsc ruutplaat
    end
    
    % for i=1:m+1   %x
    %     for j=1:n+1   %y
    %         if i==n/2 & j==m/2, textm(i,j)
    %             
    
    %   Võrrandisüsteemi koostamine --- operaator B
    va=textm(3,3);  %võrrandite arv
    ma=textm(1,3);  %muutujate arv algses süsteemis
    vs=zeros(va,ma);
    i=N/2+1;  j=M/2+1; % Keskmine sõlm
    %  ========================VS
    iv = textm(i,j);  % Võrrandi järjek. nr.
    jt = iv; 
    vs(iv,jt)=vs(iv,jt)+ 6 + 8*al^2 + 6*al^4;
    ii=i;   
    for jj=[j-2 j+2]
        jt = textm(ii,jj);  vs(iv,jt)=vs(iv,jt)+ al^4;
    end
    for jj=[j-1 j+1]
        jt = textm(ii,jj);  vs(iv,jt)=vs(iv,jt)-4*al^2*(1+al^2);
    end
    jj=j;
    for ii=[i-1 i+1]
        jt = textm(ii,jj);  vs(iv,jt)=vs(iv,jt)-4*al^2*(1+al^2);
    end
    jj=j;
    for ii=[i-2 i+2]
        jt = textm(ii,jj);  vs(iv,jt)=vs(iv,jt)+1;
    end
    for jj=[j-1 j+1]
        for ii=[i-1,i+1]
            jt = textm(ii,jj);  vs(iv,jt)=vs(iv,jt)+ 2*al^2;
        end
    end
    %  ========================VS
    
    NN=1;
    for i=M/2 : -1 : 3
        for j = M/2+1:-1:M/2-NN+1
            %  ========================VS
            iv = textm(i,j);  % Võrrandi järjek. nr.
            jt = iv; 
            vs(iv,jt)=vs(iv,jt)+ 6 + 8*al^2 + 6*al^4;
            ii=i;   
            for jj=[j-2 j+2]
                jt = textm(ii,jj);  vs(iv,jt)=vs(iv,jt)+ al^4;
            end
            for jj=[j-1 j+1]
                jt = textm(ii,jj);  vs(iv,jt)=vs(iv,jt)-4*al^2*(1+al^2);
            end
            jj=j;
            for ii=[i-1 i+1]
                jt = textm(ii,jj);  vs(iv,jt)=vs(iv,jt)-4*al^2*(1+al^2);
            end
            jj=j;
            for ii=[i-2 i+2]
                jt = textm(ii,jj);  vs(iv,jt)=vs(iv,jt)+1;
            end
            for jj=[j-1 j+1]
                for ii=[i-1,i+1]
                    jt = textm(ii,jj);  vs(iv,jt)=vs(iv,jt)+ 2*al^2;
                end
            end
            %  ========================VS
        end
        NN=NN+1;
    end
%         if kgr==1,    vs,   end
    
    % Rajatingimused
    if rajat==1,
        % 1. Jäik kinnitus:
        k1= 0.0138; %tbl. väärtus
        tittxt='Jäik kinnitus';
        for j=M/2+1:-1:3
            jt=textm(3,j);   jt1=textm(1,j);
            vs(:,jt)=vs(:,jt)+vs(:,jt1);
        end
    elseif rajat==2,
        % 2. Vaba toetus:
        k1=0.0443; %tbl. väärtus
        tittxt='Vaba toetus';
        for j=M/2+1:-1:3
            jt=textm(3,j);   jt1=textm(1,j);
            vs(:,jt)=vs(:,jt)-vs(:,jt1);
        end
    else,'rajating. pole fix',
    end
    vs=vs(:,1:va);
            if kgr==2,    vs,   end
    W=inv(vs)*ones(va,1);
            if kgr==2,    W,   end

    nu=0.3; ktbl=n^4/12/(1-nu^2);
    w=W/ktbl;
    k1w=w(1)
    mv=[mv m];
    k1wv=[k1wv k1w];
    M1= (2.6*W(1)-2.6*W(2)),    M1=M1/n^2
    if rajat==2, M4=0,
    else
        M4= -2*W(2), M4=M4/n^2,      % n=4
%         M4= -2*W(16), M4=M4/n^2,   % n=12
%         M4= -2*W(67), M4=M4/n^2,   % n=24
    end
end %m=
figure(2)
plot(mv,k1wv,mv,0*mv+k1,'r','linewidth',2),
xlabel('n'),ylabel('k_1^* ja k_1') ,title(tittxt)
times12
grid
if rajat==1, 
    print -depsc jaik_n_k_viga
end
if rajat==2, 
    print -depsc vabat_n_k_viga
end

figure(3)
viga_R=abs(k1wv-k1)/k1*100;
plot(mv,viga_R,'linewidth',2),
xlabel('n'),ylabel('k_1^* relatiivne viga, %') ,title(tittxt)
times12
grid
if rajat==2, 
    print -depsc vabat_n_k_rviga
end
if rajat==1, 
    print -depsc jaik_n_k_rviga
end

