##Naide2_7.m
#  vt lk 67
close all 
clear
## 
##dbwhere(860)
#
##=========================================================================
## This Program is written by Andres Lahe,   2015-05-10
##                    e-mail: andres.lahe@ttu.ee
## LAST MODIFIED: Andres Lahe,    2015-05-12
## Copyright (c)  2015  by Tallinn University of Technology 
##                Department of Mechanics
##                http://www.ttu.ee/%
##
## This program is free software; you can redistribute it and/or
## modify it under the terms of the GNU General Public License
## as published by the Free Software Foundation; either version 2
## of the License, or (at your option) any later version.

## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
## GNU General Public License for more details.

## You should have received a copy of the GNU General Public License
## along with this program; if not, write to the Free Software
## Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
## http://www.gnu.org/licenses/old-licenses/gpl-2.0.html
##=========================================================================
#

 global spA
 global X
 
NNK=8 
#--------- Define a matrix spA to be sparse -----------
disp('-------- Sparse matrix instantiation --------  ')
disp('  spA=sparse(NNK,NNK)  ')
spA=sparse(NNK,NNK) # nullistan võrrandisüsteemi kordajate maatriksi
#--------- Define a vector B ------------
disp('-------- Right-hand side of the equations (RHS). --------  ')
disp('  Bvb=zeros(NNK,1);  ')
Bvb=zeros(NNK,1);   # nullistan võrrandisüsteemi vabaliikmed
disp('   ')
#######
#baasi0=1.0   # baasjäikus
l=400   # cm   # varda pikkus     
#x=0.0  # cm   # koordinaat x
 # N a=0.0   # cm   # pikijõud algus
 # N mx=-150 # N    # pikijõud 
 # omega=30  # cm2
By=3.0e+06; # Ncm2 
#
GIt=1.1874*10^10 # Ncm2      # väändejäikus
#baasi0=GIt   # baasjäikus
baasi0=1.0*10^10   # baasjäikus
EIw = 2.9686*10^15 #Ncm4 # kooldumisjäikus
#######
spvF=yspWGvfmhvI(baasi0,l,GIt,EIw);  # hõreda laiendatud ülekandemaatriksi arvutus
##vB=yzWGmx(baasi0,l,l,a,mx,GIt,EIw); # koormusvektori arvutus
IIv=1;
IJv=1;
#vB=vB.*(-1);
spA=spInsertBtoA(spA,IIv,IJv,spvF);     # astab ülekandemaatriksi võrrandisüsteemi
#Bvb=InsertBtoA(Bvb,NNK,1,IIv,1,vB,4,1); # astab koormusvektori vabaliikmesse
####### Rajatingimused 
spA=spSisestaArv(spA,5,1,1);  # $theta_A$ - väändenurk
spA=spSisestaArv(spA,6,3,1);  # $B_{L}$ -bimoment 
Bvb(6,1)=-By; ## -3.0e+06;
spA=spSisestaArv(spA,7,5,1);  # $theta_A$ - väändenurk
spA=spSisestaArv(spA,8,7,1);  # $B_{L}$ -bimoment 
Bvb(8,1)= By;  ## 3.0e+06;
#  vastavd  vabaliikmed Bvb on juba nullistatud
#
spA_rank = sprank(spA) # võrrandisüsteemi astak
#A=full(spA)
spA  # hõre laiendatud ¨ulekandemaatriks
Bvb
# Teine märgikokkulepe.Väändnurk skaleeritud
Z=spA\Bvb; # lahendab võrrandisõteemist spA*X=Bvb
Z
disp('skaleerimata väändenurk skaleerimata väändenurk varda lõpus ')
theta=Z(5,1)/baasi0  # skaleerimata väändenurk varda lõpus

NEARV=1;

AlgPar(1:4,1)=[Z(1,1)/baasi0  Z(2,1) Z(3,1)  Z(4,1)];

AP=AlgPar(:,1)
baasi0=1.0
Nmitmeks=4
xx=0;
xsamm=l/Nmitmeks;
for ij=1:Nmitmeks+1  # 5 - displacements and forces at x=0.0   
Xloikes(ij,1)=xx;
   vvF=ylWGfhlin(baasi0,l,xx,GIt,EIw);

 Fvv(1:4,ij)=vvF*AP ; # +vvB;
 Fvv(5,ij)=Fvv(2,ij)+Fvv(4,ij);
xx=xx+xsamm;
endfor 
#Fvv(1,:)
Fvv;
kt=sqrt(GIt/EIw);
k=kt
  nu=l*sqrt(GIt/EIw);
 nul=nu/l;
 lnu=l/nu;
# nux=nul*x;
  nux=nul*l;
 ch=cosh(nux);
 sh=sinh(nux);
#
  
  suurused=[' theta -'; 
            '    Tt -'; 
            '     B -'; 
            '    Tw -'; 
            '  Tsum -']; 
#disp(' size(suurused) ')            
  size(suurused);          

#for i=1:5
disp(sprintf('%9s %9.2f    %9.2f     %9.2f    %9.2f    %9.2f', ' x= ',  Xloikes(1,1), Xloikes(2,1), Xloikes(3,1), Xloikes(4,1),Xloikes(5,1)))
#endfor
 #  
#
for i=1:5 
disp(sprintf('%8s %12.3e %12.3e %12.3e %12.3e %12.3e',suurused(i,:), Fvv(i,1), Fvv(i,2), Fvv(i,3), Fvv(i,4), Fvv(i,5)))
endfor
 # 
disp('  ') 
#%%%%%%%%%%%%%%%%%%%%%%%%%%%
#figure(1)
#spy(spA)
#title('spy(spA) - the sparse matrix spA(8,8) non zero elements [31%] ')
#refresh  
##print('-dfig','-landscape','spESTNaide2_1_sparse_matrix.fig');  
#refresh
#disp('-----------------------------------')
#disp('  ')
disp(' The END ')
#%%%%%%%%%%%%%%%  
  #  
 MiMaSg=minMaxGrfVnurk(Fvv);
disp('  ') 
Nmitmeks=800
xx=0;
xsamm=l/Nmitmeks;
for ij=1:Nmitmeks+1  # 5 - displacements and forces at x=0.0   
Xloikes(ij,1)=xx;
   vvF=ylWGfhlin(baasi0,l,xx,GIt,EIw);
# vvF=ylTfhlin(1.0,xx,GAr,EI);               # ylfhlin(1.0,xx,EA,GAr,EI);  
%%  vvB=yzWGmx(baasi0,l,xx,a,mx,GIt,EIw); 
#  vvB=ESTtalaKrmus(1.0,xx,Li,Fjoud,qkoormus,EI);
  
 Fvv(1:4,ij)=vvF*AP;  %% +vvB;
 Fvv(5,ij)=Fvv(2,ij)+Fvv(4,ij);
xx=xx+xsamm;
endfor 
#%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(1)
spy(spA)
title('spy(spA) - the sparse matrix spA(8,8), 20 non zero elements [31%] ')
refresh  
#print('-dfig','-landscape','spESTNaide2_1_sparse_matrix.fig');  
refresh
disp('-----------------------------------')
disp('  ')
#disp(' The END ')
#%%%%%%%%%%%%%%%  
% 
disp('  ') 
disp('  ') 
%%%%%%%%%%%  
  MiMaTa=minMaxGrfVnurk(Fvv);
  	maxthetagr=MiMaTa(1,1);
	minthetagr=MiMaTa(1,2);
	maxTtgr=MiMaTa(2,1);
	minTtgr=MiMaTa(2,2);
	maxBwgr=MiMaTa(3,1);
	minBwgr=MiMaTa(3,2);
	maxTwgr=MiMaTa(4,1);
	minTwgr=MiMaTa(4,2);        
	maxTsumgr=MiMaTa(5,1);
	minTsumgr=MiMaTa(5,2); 
%%%%%%%%%
L=l;
JA=-0.1*L;
JL=l+0.1*L;
jaotT=L/8;
XL1=L;
Xg=[0.0; XL1];
Yg=[0.0; 0.0];
%%%%%%%%%
XloikesT= Xloikes';
 FvvT=Fvv; 
%%
 XA0=XloikesT(1,1);
 YA0=0.0;
%% Xloikes(ij,1)%
 ViimaneXm=size(Xloikes);
 ViimaneX=ViimaneXm(1,1);
%%%%%%%%% 
 YAY=FvvT(1,1);
 YLY=FvvT(1,ViimaneX);
 X0=[0.0; 0.0];
 Y0=[0.0;  YAY];
 YAY=FvvT(1,1);
 XL=[L; L];
 YL=[0.0; YLY];
%%%%%%%%%   
figure(2)
%
 plot( X0,Y0,"3",XloikesT(1,:),FvvT(1,:),"3",Xg,Yg,"3",XL,YL,"3")
grid("on")
% http://www.network-theory.co.uk/docs/octave3/octave_160.html
%grid("off")
%
ax = gca ();
set (ax, "outerposition", [ 0  0  1 0.55])
%minthetagr-150maxthetagr+150
 axis([JA JL minthetagr-0.000008 maxthetagr+0.000025], "ij")  % telgede piirid [Xalg Xlõpp Zalg Zlõpp]
 ax=gca();
   xticks(1,:)=0.0:jaotT:L;
 set (ax, "xtick", xticks)
title('V66ndenurk theta ',"fontsize",10)
        xlabel('x [m]',"fontsize",10)
	ylabel('theta [rad]',"fontsize",10)
%
%%print('-deps','V66ndenurk_theta.eps');
%%print('-dpng','V66ndenurk_theta.png');
%print('-dfig','-landscape','-mono','-solid','V66ndenurk_theta.fig');
%%print('-dfig','-solid','-landscape','-color','V66ndenurk_theta.fig');
%%http://www.gnu.org/software/octave/doc/interpreter/Printing-Plots.html
refresh
%title('Kooldumisv66ndemoment Tw ',"fontsize",10)
%Vabav66ndemoment
%Koguv66ndemoment
%%%%%%%%%    
%%%%%%%%% 
 YAY=FvvT(2,1);
 YLY=FvvT(2,ViimaneX);
 X0=[0.0; 0.0];
 Y0=[0.0;  YAY];
 YAY=FvvT(2,1);
 XL=[L; L];
 YL=[0.0; YLY];
%%%%%%%%%  
figure(3)
%
 plot( X0,Y0,"3",XloikesT(1,:),FvvT(2,:),"3",Xg,Yg,"3",XL,YL,"3")
grid("on")
% http://www.network-theory.co.uk/docs/octave3/octave_160.html
%grid("off")
%
ax = gca ();
set (ax, "outerposition", [ 0  0  1 0.55])
%
 axis([JA JL minTtgr-1000 maxTtgr+1000], "ij")  % telgede piirid [Xalg Xlõpp Zalg Zlõpp]
 ax=gca();
   xticks(1,:)=0.0:jaotT:L;
 set (ax, "xtick", xticks)
title('Vabav66ndemoment Tt ',"fontsize",10)
        xlabel('x [m]',"fontsize",10)
	ylabel('Tt [Ncm]',"fontsize",10)
%
%%print('-deps','Vabav66ndemoment_Tt.eps');
%%print('-dpng','Vabav66ndemoment_Tt.png');
%print('-dfig','-landscape','-mono','-solid','Vabav66ndemoment_Tt.fig');
%%print('-dfig','-solid','-landscape','-color','Vabav66ndemoment_Tt.fig');
%%http://www.gnu.org/software/octave/doc/interpreter/Printing-Plots.html
refresh
%title('Kooldumisv66ndemoment Tw ',"fontsize",10)
%Vabav66ndemoment

%%%%%%%%% 
 YAY=FvvT(3,1);
 YLY=FvvT(3,ViimaneX);
 X0=[0.0; 0.0];
 Y0=[0.0;  YAY];
 XL=[L; L];
 YL=[0.0; YLY];
%%%%%%%%% 
figure(4)
%
 plot( X0,Y0,"3",XloikesT(1,:),FvvT(3,:),"3",Xg,Yg,"3",XL,YL,"3")
grid("on")
% http://www.network-theory.co.uk/docs/octave3/octave_160.html
%grid("off")
%
ax = gca ();
%set (ax, "outerposition", [ 0  0  1 0.4])
set (ax, "outerposition", [ 0  0  1 0.55])
%minBwgr-3800000 maxBwgr+1000000
 axis([JA JL minBwgr-3800000 maxBwgr+1000000], "ij")  % telgede piirid [Xalg Xlõpp Zalg Zlõpp]
 ax=gca();
   xticks(1,:)=0.0:jaotT:L;
 set (ax, "xtick", xticks)
title('Bimoment Bw ',"fontsize",10)
        xlabel('x [m]',"fontsize",10)
	ylabel('B [Ncm2]',"fontsize",10)
%
%%print('-deps','Bimoment_Bw.eps');
%%print('-dpng','Bimoment_Bw.png');
%print('-dfig','-landscape','-mono','-solid','Bimoment_Bw.fig');
%%print('-dfig','-solid','-landscape','-color','Bimoment_Bw.fig');
%%http://www.gnu.org/software/octave/doc/interpreter/Printing-Plots.html
refresh
% 
%%%%%%%%% 
 YAY=FvvT(4,1);
 YLY=FvvT(4,ViimaneX);
 X0=[0.0; 0.0];
 Y0=[0.0;  YAY];
 YAY=FvvT(4,1);
 XL=[L; L];
 YL=[0.0; YLY];
%%%%%%%%%  
figure(5)
%
 plot( X0,Y0,"3",XloikesT(1,:),FvvT(4,:),"3",Xg,Yg,"3",XL,YL,"3")
grid("on")
% http://www.network-theory.co.uk/docs/octave3/octave_160.html
%grid("off")
%
ax = gca ();
set (ax, "outerposition", [ 0  0  1 0.55])
%minTwgr-4000 maxTwgr+1500
 axis([JA JL minTwgr-4000 maxTwgr+1500], "ij")  % telgede piirid [Xalg Xlõpp Zalg Zlõpp]
 ax=gca();
   xticks(1,:)=0.0:jaotT:L;
 set (ax, "xtick", xticks)
title('Kooldumisv66ndemoment Tw ',"fontsize",10)
        xlabel('x [m]',"fontsize",10)
	ylabel('Tw [Ncm]',"fontsize",10)
%
%%print('-deps','Kooldumisv66ndemoment_Tw.eps');
%%print('-dpng','Kooldumisv66ndemoment_Tw.png');
%print('-dfig','-landscape','-mono','-solid','Kooldumisv66ndemoment_Tw.fig');
%%print('-dfig','-solid','-landscape','-color','Kooldumisv66ndemoment_Tw.fig');
%%http://www.gnu.org/software/octave/doc/interpreter/Printing-Plots.html
refresh
%title('Kooldumisv66ndemoment Tw ',"fontsize",10)
%Vabav66ndemoment

%%%%%%%%% 
 YAY=FvvT(5,1);
 YLY=FvvT(5,ViimaneX);
 X0=[0.0; 0.0];
 Y0=[0.0;  YAY];
 YAY=FvvT(5,1);
 XL=[L; L];
 YL=[0.0; YLY];
%%%%%%%%%   
figure(6)
%
 plot( X0,Y0,"3",XloikesT(1,:),FvvT(5,:),"3",Xg,Yg,"3",XL,YL,"3")
grid("on")
% http://www.network-theory.co.uk/docs/octave3/octave_160.html
%grid("off")
%
ax = gca ();
set (ax, "outerposition", [ 0  0  1 0.55])
%minTwgr-2000.0 maxTwgr+100.0
 axis([JA JL -0.000002 0.000002], "ij")  % telgede piirid [Xalg Xlõpp Zalg Zlõpp]
 ax=gca();
   xticks(1,:)=0.0:jaotT:L;
 set (ax, "xtick", xticks)
title('Koguv66ndemoment Tsum ',"fontsize",10)
        xlabel('x [m]',"fontsize",10)
	ylabel('Tsum [Ncm]',"fontsize",10)
%
%%print('-deps','Koguv66ndemoment_Tsum.eps');
%%print('-dpng','Koguv66ndemoment_Tsum.png');
%print('-dfig','-landscape','-mono','-solid','Koguv66ndemoment_Tsum.fig');
%%print('-dfig','-solid','-landscape','-color','Koguv66ndemoment_Tsum.fig');
%%http://www.gnu.org/software/octave/doc/interpreter/Printing-Plots.html
refresh
%title('Kooldumisv66ndemoment Tw ',"fontsize",10)
%Vabav66ndemoment
%Koguv66ndemoment
%%%%%%%%% 
    
  
  
  
  
  
  
  
   
   
   
   
   
   
   
   
   
   
   
   








