%......沒用迭代計算擾度的程序..........%
%......D為變值,g為改進過..........%
function Calculate_g(row,column,long,wide,E,h,p,coefficient,force,a,EI1,EI2,EI3,EI4)
% 輸入參數，row,column為橫縱方向上的節點數，長，寬，彈性模量，板厚，泊松比，地基系數,集中力
h1=long/(row-1);%X方向上的步長
h2=wide/(column-1);%Y方向上的步長
t1=E*h^3/(12*(1-p^2));%暫存值
D=zeros(row,column);%賦值
for i=1:row
    for j=1:column
        D(i,j)=t1*(1-4*(1-a)*(((i-1)*h1/long-1/2)^2+((j-1)*h2/wide-1/2)^2));
    end
end
DL=zeros(row-1,column-1);%賦值
ul=zeros(row-1,column-1);%賦值
for i=1:(row-1)
    for j=1:(column-1)
        DL(i,j)=t1*(1-4*(1-a)*(((i-0.5)*h1/long-1/2)^2+((j-0.5)*h2/wide-1/2)^2));
        ul(i,j)=p;
    end
end
u=zeros(row,column);
C=zeros(row,column);
for i=1:row
    for j=1:column
        C(i,j)=coefficient;%地基系數
        u(i,j)=p;
    end
end
t2=force/(h1*h2);%計算作用力
count=row*column;%count為方程的系數矩陣的維數，行*列%
%...................調用擾度函數.................................. %
DiJi(row,column,h1,h2,D,DL,u,ul,t2,C,EI1,EI2,EI3,EI4,count)
%梁1,2,3,4的取值，初值設為0，常數 %
%........................................................
%計算擾度的函數,A為系數矩陣,常數項b,X為擾度,row行節點數,column為列節點數,h1,h2為步長,u為常數
%D,u，c分別為各點的彎曲剛度值與...，需要輸入size(D)=row*column維矩陣，u，c維數也一樣
%DL為節點之中間的彎曲剛度值，輸入size(DL)=row*column，ul也一樣
%q為輸入值，size(q)=row*column,1%
%                row=N1+1，column=N2+1          %
function DiJi(row,column,h1,h2,D,DL,u,ul,t2,C,EI1,EI2,EI3,EI4,count)
%D代表剛度系數，DL代表/D的剛度系數%
%    計算方程的系數矩陣A，標號從(1,1)算起  %
bc0=4;%A對角線上的步長%
bc1=1;%A對角線上一行的步長%
bc2=1;%A對角線上兩行的步長%
bc3=4;%S對角線的步長%
bc4=1;%S對角線上一行的步長%
bc5=1;%S對角線下一行的步長%
bc6=4;%T對角線的步長%
h12=h1^2*h2^2;
k1=(h2/h1)^2;
k2=1/k1;
N1=row-1;%N1為地基最右端的坐標,N1>=7，行坐標
N2=column-1;%N2為地基最上端的坐標,N2>=7，列坐標
q=zeros(count,1);%常數項                                                                                                      
b=zeros(count,1);%常數項
%.........地基為N1*N2 ....需改動A,B,C,D點的位置...
q(column*ceil(N1/4)+ceil(column/4),1)=t2;%A點作用力,X1，X2坐標為長(N1),高(N2)的的四分子一，
q(column*ceil(N1*3/4)+ceil(column/4),1)=t2;%B點作用力，X1坐標為長(N1),高(N2)的四分子三
q(column*ceil(N1/4)+ceil(column*3/4),1)=t2;%C點作用力，X2坐標為ceil(N1/4)=5
q(column*ceil(N1*3/4)+ceil(column*3/4),1)=t2;%D點作用力，X1，X2坐標為長(N1),高(N2)的的四分子三
%.............對q,b稀疏化，節省內存.........................
q=sparse(q);
b=sparse(b);
%.......................................................
A=zeros(count);%生成n*n維的零矩陣，以便對其賦值%
%.............對A稀疏化，節省內存........................
A=sparse(A);
%......................................................
%原文中的DL(0.5,0.5)為DL(1,1)勿需改動,D(1,0)為D(2,1),u(0,1)為u(1,2),ul與原文相同，C(0,0)為C(1,1)，循環時DL，ul跟D，C都得變化,但是x坐標不變，y坐標變時不加1，D，C的x，y坐標都在原來基礎上加1%
%                 輸入A中的數據元素：對角線上一行的元素         %
%             下面為地基中的第一列               %
  A(1,2)=-2*DL(1,1)*(1-ul(1,1))-k2*(D(1,2)*(1-u(1,2)*u(1,2))+2*EI1/h1)-5*EI1*k2/h1;
%a'(0,1)=-2*DL(1,1)*(1-ul(1,1))-k2*(D(0,1)*(1-u(0,1)*u(0,1))+2*EI1/h1)-5*EI1*k2/h1%
%輸入紙1，/a01,W01的系數%
%點(0,1)的D(0,1)相當于矩陣中的D(1,2)%
  A(2,3)=-2*DL(1,2)*(1-ul(1,2))-k2*(D(1,2)*(1-u(1,2)*u(1,2))+2*EI1/h1)-k2*(D(1,3)*(1-u(1,3)*u(1,3))+2*EI1/h1)-10*EI1*k2/h1;
%a'(0,2)=-2*DL(1,2)*(1-ul(1,2))-k2*(D(0,1)*(1-u(0,1)*u(0,1))+2*EI1/h1)-k2*(D(0,2)*(1-u(0,2)*u(0,2))+2*EI1/h1)-10*EI1*k2/h1%
%輸入紙2，/a02,W02的系數%
  A(3,4)=-2*DL(1,3)*(1-ul(1,3))-k2*(D(1,3)*(1-u(1,3)*u(1,3))+2*EI1/h1)-k2*(D(1,4)*(1-u(1,4)*u(1,4))+2*EI1/h1)-2*EI1*k2/h1;
%a'(0,3)=-2*DL(1,3)*(1-ul(1,3))-k2*(D(0,2)*(1-u(0,2)*u(0,2))+2*EI1/h1)-k2*(D(0,3)*(1-u(0,3)*u(0,3))+2*EI1/h1)-2*EI1*k2/h1%
%輸入紙3，/a03,W03的系數%
  A(4,5)=-2*DL(1,4)*(1-ul(1,4))-k2*(D(1,4)*(1-u(1,4)*u(1,4))+2*EI1/h1)-k2*(D(1,5)*(1-u(1,5)*u(1,5))+2*EI1/h1);
%a'(0,4)=-2*DL(1,4)*(1-ul(1,4))-k2*(D(0,3)*(1-u(0,3)*u(0,3))+2*EI1/h1)-k2*(D(0,4)*(1-u(0,4)*u(0,4))+2*EI1/h1)%
%輸入紙4，/a04,W04的系數%
%N2為列節點末點的坐標，通式%
for j=5:N2-3  
    A(j,j+1)=-2*DL(1,j)*(1-ul(1,j))-k2*(D(1,j)*(1-u(1,j)*u(1,j))+2*EI1/h1)-k2*(D(1,j+1)*(1-u(1,j+1)*u(1,j+1))+2*EI1/h1);
    %a'(0,j)=-2*DL(1,j)*(1-ul(1,j))-k2*(D(0,j-1)*(1-u(0,j-1)*u(0,j-1))+2*EI1/h1)-k2*(D(0,j)*(1-u(0,j)*u(0,j))+2*EI1/h1)%
    %輸入紙5，/a0k,W0(j+1)的系數,其中j=5,6,...,N2-3%
end
A(N2-2,N2-1)=-2*DL(1,N2-2)*(1-ul(1,N2-2))-k2*(D(1,N2-2)*(1-u(1,N2-2)*u(1,N2-2))+2*EI1/h1)-k2*(D(1,N2-1)*(1-u(1,N2-1)*u(1,N2-1))+2*EI1/h1)-2*EI1*k2/h1;
% a'(0,N2-2)=-2*DL(1,N2-2)*(1-ul(1,N2-2))-k2*(D(0,N2-3)*(1-u(0,N2-3)*u(0,N2-3))+2*EI1/h1)-k2*(D(0,N2-2)*(1-u(0,N2-2)*u(0,N2-2))+2*EI1/h1)-2*EI1*k2/h1%
%輸入紙6，/a0(N2-2),W0(N2-2)的系數%
 A(N2-1,N2)=-2*DL(1,N2-1)*(1-ul(1,N2-1))-k2*(D(1,N2-1)*(1-u(1,N2-1)*u(1,N2-1))+2*EI1/h1)-k2*(D(1,N2)*(1-u(1,N2)*u(1,N2))+2*EI1/h1)-10*EI1*k2/h1;
%a'(0,N2-1)=-2*DL(1,N2-1)*(1-ul(1,N2-1))-k2*(D(0,N2-2)*(1-u(0,N2-2)*u(0,N2-2))+2*EI1/h1)-k2*(D(0,N2-1)*(1-u(0,N2-1)*u(0,N2-1))+2*EI1/h1)-10*EI1*k2/h1%
%輸入紙7，/a0(N2-1),W0(N2-1)的系數%
A(N2,column)=-2*DL(1,N2)*(1-ul(1,N2))-k2*(D(1,N2)*(1-u(1,N2)*u(1,N2))+2*EI1/h1)-5*EI1*k2/h1;
% a'(0,N2)=-2*DL(1,N2)*(1-ul(1,N2))-k2*(D(0,N2-1)*(1-u(0,N2-1)*u(0,N2-1))+2*EI1/h1)-5*EI1*k2/h1%
%輸入紙8，/a0N2,W0N2的系數%

%             下面為地基中的第二列到N1列 i=1,2...N1-1               %
for i=N2+2:column:((N1-1)*column+1)          %從N2+2處開始,步長為N2+1%
    A(i,i+1)=-2*D(i-bc1*N2,2)*(k2+u(i-bc1*N2,2))-2*DL(i-bc1*N2-1,1)*(1-ul(i-bc1*N2-1,1))-2*DL(i-bc1*N2,1)*(1-ul(i-bc1*N2,1));
    %a'(i,1)=-2*D(i,1)*(k2+u(i,1))              -2*DL(i,1)*(1-ul(i,1))-2*DL(i+1,1)*(1-ul(i+1,1))%
    %輸入論文55，/ai.1%
    for j=i+1:i+N2-2
        A(j,j+1)=-2*D(i-bc1*N2,j-i+1)*(k2+u(i-bc1*N2,j-i+1))-2*D(i-bc1*N2,j-i+2)*(k2+u(i-bc1*N2,j-i+2))-2*DL(i-bc1*N2-1,j-i+1)*(1-ul(i-bc1*N2-1,j-i+1))-2*DL(i-bc1*N2,j-i+1)*(1-ul(i-bc1*N2,j-i+1));
      %  a'(i,j)=-2*D(i,j-1)*(k2+u(i,j-1))                  -2*D(i,j)*(k2+u(i,j))                      -2*DL(i,j)*(1-ul(i,j))                          -2*DL(i+1,j)*(1-ul(i+1,j))%
        %輸入論文55,/aij，其中j=2,...,N2-1%
    end
    A(i+N2-1,i+N2)=-2*D(i-bc1*N2,N2)*(k2+u(i-bc1*N2,N2))-2*DL(i-bc1*N2-1,N2)*(1-ul(i-bc1*N2-1,N2))-2*DL(i-bc1*N2,N2)*(1-ul(i-bc1*N2,N2));
   %      a'(i,N2)=-2*D(i,N2-1)*(k2+u(i,N2-1))        -2*DL(i,N2)*(1-ul(i,N2))                -2*DL(i+1,N2)*(1-ul(i+1,N2))%
    %輸入論文55,/aiN2%
    bc1=bc1+1;
end

%             下面為地基中的第N1+1列               %
  A(N1*column+1,N1*column+2)=-2*DL(N1,1)*(1-ul(N1,1))-k2*(D(row,2)*(1-u(row,2)*u(row,2))+2*EI3/h1)-5*EI3*k2/h1;
%                   a'(N1,1)=-2*DL(N1,1)*(1-ul(N1,1))-k2*(D(N1,1)*(1-u(N1,1)*u(N1,1))+2*EI3/h1)-5*EI3*k2/h1%
%輸入紙24，/aN1.1,WN1.1的系數%
A(N1*column+2,N1*column+3)=-2*DL(N1,2)*(1-ul(N1,2))-k2*(D(row,2)*(1-u(row,2)*u(row,2))+2*EI3/h1)-k2*(D(row,3)*(1-u(row,3)*u(row,3))+2*EI3/h1)-10*EI3*k2/h1;
%                 a'(N1,2)=-2*DL(N1,2)*(1-ul(N1,2))-k2*(D(N1,1)*(1-u(N1,1)*u(N1,1))+2*EI3/h1)      -k2*(D(N1,2)*(1-u(N1,2)*u(N1,2))+2*EI3/h1)-10*EI3*k2/h1%
%輸入紙25，/aN1.2,WN1.2的系數%
A(N1*column+3,N1*column+4)=-2*DL(N1,3)*(1-ul(N1,3))-k2*(D(row,3)*(1-u(row,3)*u(row,3))+2*EI3/h1)-k2*(D(row,4)*(1-u(row,4)*u(row,4))+2*EI3/h1)-2*EI3*k2/h1;
%                 a'(N1,3)=-2*DL(N1,3)*(1-ul(N1,3))-k2*(D(N1,2)*(1-u(N1,2)*u(N1,2))+2*EI3/h1)-k2*(D(N1,3)*(1-u(N1,3)*u(N1,3))+2*EI3/h1)-2*EI3*k2/h1%
%輸入紙26，/aN1.3,WN1.3的系數%
A(N1*column+4,N1*column+5)=-2*DL(N1,4)*(1-ul(N1,4))-k2*(D(row,4)*(1-u(row,4)*u(row,4))+2*EI3/h1)-k2*(D(row,5)*(1-u(row,5)*u(row,5))+2*EI3/h1);
%                 a'(N1,4)=-2*DL(N1,4)*(1-ul(N1,4))-k2*(D(N1,3)*(1-u(N1,3)*u(N1,3))+2*EI3/h1)-k2*(D(N1,4)*(1-u(N1,4)*u(N1,4))+2*EI3/h1)%
%輸入紙27，/aN1.4,WN1.4的系數%
for j=N1*column+5:row*column-4
    A(j,j+1)=-2*DL(N1,j-N1*column)*(1-ul(N1,j-N1*column))-k2*(D(row,j-N1*column)*(1-u(row,j-N1*column)*u(row,j-N1*column))+2*EI3/h1)-k2*(D(row,j-N1*column+1)*(1-u(row,j-N1*column+1)*u(row,j-N1*column+1))+2*EI3/h1);
    %a'(N1,j)=-2*DL(N1,j)*(1-ul(N1,j))                   -k2*(D(N1,j-1)*(1-u(N1,j-1)*u(N1,j-1))+2*EI3/h1)                              -k2*(D(N1,j)*(1-u(N1,j)*u(N1,j))+2*EI3/h1)%
    %輸入紙28,/aN1.k，WN1.(j+1)的系數,其中K=5,6,...,N2-3%
end
A(row*column-3,row*column-2)=-2*DL(N1,N2-2)*(1-ul(N1,N2-2))-k2*(D(row,N2-2)*(1-u(row,N2-2)*u(row,N2-2))+2*EI3/h1)-k2*(D(row,N2-1)*(1-u(row,N2-1)*u(row,N2-1))+2*EI3/h1)-2*EI3*k2/h1;
%                      a'(N1,N2-2)=-2*DL(N1,N2-2)*(1-ul(N1,N2-2))-k2*(D(N1,N2-3)*(1-u(N1,N2-3)*u(N1,N2-3))+2*EI3/h1)      -k2*(D(N1,N2-2)*(1-u(N1,N2-2)*u(N1,N2-2))+2*EI3/h1)-2*EI3*k2/h1%
%輸入紙29，/aN1.(N2-2),WN1.(N2-2)的系數%
A(row*column-2,row*column-1)=-2*DL(N1,N2-1)*(1-ul(N1,N2-1))-k2*(D(row,N2-1)*(1-u(row,N2-1)*u(row,N2-1))+2*EI3/h1)-k2*(D(row,N2)*(1-u(row,N2)*u(row,N2))+2*EI3/h1)-10*EI3*k2/h1;
%                      a'(N1,N2-1)=-2*DL(N1,N2-1)*(1-ul(N1,N2-1))-k2*(D(N1,N2-2)*(1-u(N1,N2-2)*u(N1,N2-2))+2*EI3/h1)      -k2*(D(N1,N2-1)*(1-u(N1,N2-1)*u(N1,N2-1))+2*EI3/h1)-10*EI3*k2/h1 %
%輸入紙30，/aN1.(N2-1),WN1.(N2-1)的系數%
A(row*column-1,row*column)=-2*DL(N1,N2)*(1-ul(N1,N2))-k2*(D(row,N2)*(1-u(row,N2)*u(row,N2))+2*EI3/h1)-5*EI3*k2/h1;
%                      a'(N1,N2)=-2*DL(N1,N2)*(1-ul(N1,N2))-k2*(D(N1,N2-1)*(1-u(N1,N2-1)*u(N1,N2-1))+2*EI3/h1)-5*EI3*k2/h1 %
%輸入紙31，/aN1.N2,WN1.N2的系數%

 %                 輸入A中的數據元素：對角線上兩行的元素         %
 %             下面為地基中的第一列               %
  A(1,3)=1/2*k2*(D(1,2)*(1-u(1,2)*u(1,2))+2*EI1/h1)+4*EI1*k2/h1;
%a"(0,1)=1/2*k2*(D(0,1)*(1-u(0,1)*u(0,1))+2*EI1/h1)+4*EI1*k2/h1%
%輸入紙1，//a0.1,W0.2的系數%
  A(2,4)=1/2*k2*(D(1,3)*(1-u(1,3)*u(1,3))+2*EI1/h1)+5/2*EI1*k2/h1;
%a"(0,2)=1/2*k2*(D(0,2)*(1-u(0,2)*u(0,2))+2*EI1/h1)+5/2*EI1*k2/h1%
%輸入紙2，//a0.2,W0.3的系數%
  A(3,5)=1/2*k2*(D(1,4)*(1-u(1,4)*u(1,4))+2*EI1/h1);
%a"(0,3)=1/2*k2*(D(0,3)*(1-u(0,3)*u(0,3))+2*EI1/h1) %