clear all

filename='blade10.txt';
fid=fopen(filename,'w+t');%this loads the file with all the data
table=(readtable('IEA_15MW_240_RWT_Blade.txt'));
X=[table.DistanceAlongZAxis];
mass=(table.Mass_unitLength);
neutrx=(table.NeutralAxis_X_);
neutry=(table.NeutralAxis_Y_);
massaxis=(table.MassAxisOrientation)*pi/180;
r=(table.RadiiOfGyrationRatio);
inertiaall=table.MassMomentOfInertia_unitLength;
thi=table.Thickness
EA=table.AxialStiffness_EA_;
E_I_XX=(table.BendingStiffnessAboutXP_EIXp_);
E_I_YY=(table.BendingStiffnessAboutYP_EIyp_);
G_J=table.TorsionalStiffness_GJ_;
G_A_xp=table.ShearStiffnessAlongXP_GAxp_;
G_A_yp=table.ShearStiffnessAlongYP_GAyp_;
beta=table.PrincipalAxisOrientation;
ch=table.Chord;
I_M=double(zeros(6,6));
I_M_1=double(zeros(6,6));
I_M_2=double(zeros(6,6));
Sx=table.ShearCentre_X__;
Sy=table.ShearCentre_Y__;
xc=table.CentreOfMass_X__;
yc=table.CentreOfMass_Y__;%starting the Mass, Inertia and Ratio interpolation
SimTime=200;%time for the simulation
tstep=0.03;%timestep
n=40
meanx=(0:117/n:117)';
m=interp1(X,mass,meanx);
inertia=interp1(X,inertiaall,meanx);
thi=interp1(X,thi,meanx)
R=interp1(X,r,meanx);
v=(0:117/n:(117))';
E_A=interp1(X,EA,v);
massaxis=interp1(X,massaxis,meanx);
E_I_XX=interp1(X,E_I_XX,v);
E_I_YY=interp1(X,E_I_YY,v);
G_J=interp1(X,G_J,v);
G_A_xp=interp1(X,G_A_xp,v);
G_A_yp=interp1(X,G_A_yp,v);
neutrx=interp1(X,neutrx,v);
neutry=interp1(X,neutry,v);
SX=interp1(X,Sx,meanx);
SY=interp1(X,Sy,meanx);
Sx=interp1(X,Sx,v)
Sy=interp1(X,Sy,v)
CH=interp1(X,ch,v)
beta=interp1(X,beta,v);
COS_beta=cosd(beta );
SIN_beta=sind(beta);
pitch=270;
mu=1
xc=interp1(X,xc,meanx);
yc=interp1(X,yc,meanx);
ch=interp1(X,ch,meanx');
fprintf(fid,'begin: data;\n');%printing the first part of the code
fprintf(fid,'problem: initial value;\n');
fprintf(fid,'end: data;\n');
fprintf(fid,'begin: initial value;\n');
fprintf(fid,'    initial time: 0.;\n');
fprintf(fid,'    final time:%d.;\n',SimTime); %total time of the simulation     
fprintf(fid,'    time step: %f;\n', tstep);%timestep
fprintf(fid,'    max iterations: 1000;\n    tolerance: %f;\nderivatives coefficient:1.e-6; derivatives tolerance:1.e7;\n output:iterations;nonlinear solver: newton raphson, true;\n',0.1 );%tolerance,derivatives data and outputs
fprintf(fid,'\nend: initial value;\n');
fprintf(fid,'begin: control data;\n ');
fprintf(fid,'default beam output:position;')
fprintf(fid,'    structural nodes:\n');%the nodes 
fprintf(fid,'		+%d     #clamped nodes\n    ',1);
fprintf(fid,'		+%d    # other nodes\n\n    ;\n',length(meanx)-1);
fprintf(fid,'    rigid bodies:\n');%the rigid bodies which includes mass and mass moments of inertia
fprintf(fid,'		+%d #mass of other nodes \n    ;\n',(length(meanx)-1)*mu);
fprintf(fid,'    joints:\n');%joints 
fprintf(fid,'		+1    # clamp \n    ;	\n');
fprintf(fid,'    beams:\n');%the beams 
fprintf(fid,'		+%d    # the whole beam \n    ;\n',length(meanx)*0.5);
fprintf(fid,'    gravity;\n forces:\n 1;');%gravity
fprintf(fid,'\nend: control data;\n');
fprintf(fid,'begin: nodes;\n');%starting the nodes analysis
fprintf(fid,'    structural: 1,static,\n');%first clamped node
fprintf(fid,'		reference, global,0,%f,%f,\n',0,0);
fprintf(fid,'		eye,\n');
fprintf(fid,'		null,\n');
fprintf(fid,'		null;\n       \n');
    fprintf(fid,'include:"blade.nod10.txt";\n');%external inclusion for nodes 
 
fprintf(fid,'end: nodes; \n begin: elements;\n');%starting the elements (beams,bodies, joints,and others) analysis
fprintf(fid,'	joint: 1, clamp, 1, node, node;\n');
fprintf(fid,'\n include:"blade.elm10.txt";');%external inclusion for bodies and beams
fprintf(fid,'\n	gravity:uniform,\n');
fprintf(fid,'0.,0.,1., 	\n');
fprintf(fid,'	const,9.81;');
fprintf(fid,'\n end: elements;');
fclose(fid);
fid=fopen('blade.nod10.txt','w+t');
K=zeros(6,6);
K_1=zeros(6,6);
K_2=zeros(6,6);
Matrix={};
for i=1:length(v) %starting to build all the  6X6 Matrices,AM is axial and SM is shear
 
    A_M(1,1)=E_A(i,1);
    A_M(1,2)=(neutry(i)-Sy(i)*ch(i)/100)*E_A(i,1);
    A_M(1,3)=-(neutrx(i)-Sx(i)*ch(i)/100)*E_A(i,1);
    S_M(1,1)=G_A_xp(i)*COS_beta(i,1)^2+G_A_yp(i)*SIN_beta(i)^2;
    S_M(1,2)=(G_A_yp(i)-G_A_xp(i))*COS_beta(i)*SIN_beta(i);
    S_M(2,2)=G_A_yp(i)*COS_beta(i,1)^2+G_A_xp(i)*SIN_beta(i)^2;
    S_M(1,3)=0;
    S_M(2,3)=0;
    S_M(3,3)=G_J(i,1);
     A_M(2,2)=(E_I_XX(i,1)*COS_beta(i,1)^2 +E_I_YY(i,1)*SIN_beta(i,1)^2)+E_A(i,1)*(neutry(i)-Sy(i)*ch(i)/100)^2;
    A_M(2,3)=(E_I_XX(i,1)-E_I_YY(i,1))*COS_beta(i,1)*SIN_beta(i,1)-(neutrx(i)-Sx(i)*ch(i)/100)*(neutry(i)-Sy(i)*ch(i)/100)*E_A(i,1);
    A_M(3,3)=E_I_XX(i,1)*SIN_beta(i,1)^2 +E_I_YY(i,1)*COS_beta(i,1)^2+E_A(i,1)*(neutrx(i)-Sx(i)*ch(i)/100)^2;
    A_M(2,1)=A_M(1,2);
    A_M(3,1)=A_M(1,3);
    A_M(3,2)=A_M(2,3);
    S_M(2,1)=S_M(1,2);
    S_M(3,1)=S_M(1,3);
    S_M(3,2)=S_M(2,3);
    K([1,5,6],[1,5,6])=A_M;
    K(2:4,2:4)=S_M;
    Matrix{i}=K;%saving all the data in a cell array
end 
for i=1:length(Matrix)
    for j=1:6
        for k=1:6
Matrix_I{i}(j,k)=Matrix{i}(j,k);
        end 
    end 
end



lengthbeam=X(end)/(length(v)-1)*2;%length of each beam
for i=1:length(meanx)
    YC(i)=yc(i)*ch(i)/100
    XC(i)=xc(i)*ch(i)/100
end 
YC=YC'
XC=XC'


for i=1:length(meanx)-1           %integrating the density (linear) in order to get the mass, required simps.m for simpson integration 
    if i<length(meanx)
iii(i)=simps(meanx(i:i+1),inertia(i:i+1))
mmm(i)=simps(meanx(i:i+1),m(i:i+1))
ZZZ(i)=simps(meanx(i:i+1),m(i:i+1).*meanx(i:i+1))/(m(i)*lengthbeam/4+m(i+1)*lengthbeam/4)
YYY(i)=simps(meanx(i:i+1),m(i:i+1).*YC(i:i+1))/mmm(i)
XXX(i)=simps(meanx(i:i+1),m(i:i+1).*XC(i:i+1))/mmm(i)
    else 
        iii(i)=simps(meanx(i:i+1),inertia(i-1:i))
        mmm(i)=simps(meanx(i:i+1),m(i-1:i))
        ZZZ(i)=simps(meanx(i:i+1),m(i:i-1).*meanx(i:i+1))/mmm(i)
        
    end 
end 
for i=1:length(ZZZ)
    static(i)=ZZZ(i)*mmm(i)
end 
ZCNEW=sum(static)/sum(mmm)

j=1;

 for j=1:length(mmm)%allocation 
     if j==1
     massmean(j)=mmm(j)
     inertiaavg(j)=iii(j)
     
     else 
         massmean(j)=mmm(j)
         inertiaavg(j)=iii(j)
     end 
 end 


 
 sum(massmean)%just to check the overall mass
   
 fi=0;
 
 %the part below is useful for the centre of mass
Area=ch.*(thi.*ch')'/100
 radius=sqrt(Area/pi)
 radiusx=(sqrt((inertiaavg./(1+R(2:end)'.^2))./massmean))'
 radiusy=(R(2:end).*radiusx)
        
for i=2:length(meanx)%printing the second part of the codes:nodes
 
fprintf(fid,'    structural: %d, dynamic, \n',i);
if i==29
fprintf(fid,'		%d,%d,%d,\n',meanx(i),0,0); 
else 
   fprintf(fid,'		%d,%d,%d,\n',meanx(i),0,0);  
end 
fprintf(fid,'		matr,1,0,0,0,%d,%d,0,%d,%d,\n',1,0,0,1);
fprintf(fid,'		null,\n'); 
fprintf(fid,'		null; \n      \n');
end
fclose(fid);
fid=fopen('blade.elm10.txt.','w+t');%printing the third part of the codes:bodies with masses, and  inertia and beams
n=1;% used for the beams later 
fprintf(fid,'force:3,absolute,%d,position,null,0,0,1,const,%d;',length(meanx),0);%external forces
q=1000;

fprintf(fid,'	body: %d , %d \n, \n',1,2);
fprintf(fid,'		+%d,\n',massmean(1));      %lenght of element times density
fprintf(fid,'    reference,node,%d,%d,%d,\n',0,yc(1)*ch(1)/100,xc(1)*ch(1)/100);         %centre of mass and below, inertia
fprintf(fid,'		diag,%f, %f, %f;\n',inertiaavg(j),massmean(1)*radiusx(1)^2,massmean(1)*radiusy(1)^2);
    c(1,1)=massmean(1)/mu
    d(1,1)=(inertiaavg(j)/(1+R(j)^2))/mu/8
    
     fprintf(fid,'inertia:%d, name,"%d", position,null,orientation,eye,body,%d,output, yes;',1,1,1);
   
     
     
     
q=q+500;
for j=2:length(massmean) %bodies start here
    
    for i=1:mu
        fprintf(fid,'	body: %d , %d \n, \n',j+q,j+1);
fprintf(fid,'		+%d,\n',massmean(j));      %lenght of element times density
fprintf(fid,'     %d,%d,%d,\n',0,yc(j+1)*ch(j+1)/100,xc(j+1)*ch(j+1)/100);         %centre of mass and below, inertia
fprintf(fid,'		diag,%f, %f, %f;\n',inertiaavg(j),radiusx(j)^2*massmean(j),radiusy(j)^2*massmean(j));
    c(i,j+1)=massmean(j)
    d(i,j+1)=(inertiaavg(j)/(1+R(j)^2))/mu+massmean(j)*(meanx(j+1)-meanx(j))^2
      fprintf(fid,'inertia:%d, name,"%d", position,null,orientation,eye,body,%d,output, yes;',j+q,j+q,j+q);
       q=q+500;
    end
   
end 
q=1000
for j=1:length(massmean)

        

q=q+500 
end 
    fprintf(fid,'inertia:%d, name,"%d", position,null,orientation,eye,body,all,output, yes;',j+q,j+q);



    
j=1;


		

            



for j=2:round(length(meanx)/2)%beams start here 
    k=2*j-2;
    
fprintf(fid,'\n	beam3: %d,\n',k/2);
fprintf(fid,"		"+(k-1)+", %d,%d,%d,    \n",0,SX(k-1)*ch(k-1)/100,SY(k-1)*ch(k-1)/100)
fprintf(fid,"		"+(k)+", %d,%d,%d,   \n",0,SX(k)*ch(k)/100,SY(k)*ch(k)/100)
fprintf(fid,"		"+(k+1)+", %d,%d,%d,    \n",0,SX(k+1)*ch(k+1)/100,SY(k+1)*ch(k+1)/100);
fprintf(fid,'		matr,1,0,0,0,%d,%d,0,%d,%d,\n',cosd(pitch),-sind(pitch),sind(pitch),cosd(pitch));
fprintf(fid,'		linear viscoelastic generic,    \n');
fprintf(fid,'			sym, %e, %e, %e, %e, %e, %e,   %e, %e, %e, %e, %e,  %e, %e, %e, %e,    %e, %e, %e,   %e,%e,  %e,\n',Matrix_I{n+1}(1,1),Matrix_I{n+1}(1,2),Matrix_I{n+1}(1,3),Matrix_I{n+1}(1,4),Matrix_I{n+1}(1,5),Matrix_I{n+1}(1,6),Matrix_I{n+1}(2,2),Matrix_I{n+1}(2,3),Matrix_I{n+1}(2,4),Matrix_I{n+1}(2,5),Matrix_I{n+1}(2,6),Matrix_I{n+1}(3,3),Matrix_I{n+1}(3,4),Matrix_I{n+1}(3,5),Matrix_I{n+1}(3,6),Matrix_I{n+1}(4,4),Matrix_I{n+1}(4,5),Matrix_I{n+1}(4,6),Matrix_I{n+1}(5,5),Matrix_I{n+1}(5,6),Matrix_I{n+1}(6,6));
fprintf(fid,'proportional,%f,',0);
fprintf(fid,'		matr,1,0,0,0,%d,%d,0,%d,%d,\n',cosd(pitch),-sind(pitch),sind(pitch),cosd(pitch))
fprintf(fid,'		linear viscoelastic generic,    \n');
fprintf(fid,'			sym, %e, %e, %e, %e, %e, %e,   %e, %e, %e, %e, %e,  %e, %e, %e, %e,    %e, %e, %e,   %e,%e,  %e,\n',Matrix_I{n+1}(1,1),Matrix_I{n+1}(1,2),Matrix_I{n+1}(1,3),Matrix_I{n+1}(1,4),Matrix_I{n+1}(1,5),Matrix_I{n+1}(1,6),Matrix_I{n+1}(2,2),Matrix_I{n+1}(2,3),Matrix_I{n+1}(2,4),Matrix_I{n+1}(2,5),Matrix_I{n+1}(2,6),Matrix_I{n+1}(3,3),Matrix_I{n+1}(3,4),Matrix_I{n+1}(3,5),Matrix_I{n+1}(3,6),Matrix_I{n+1}(4,4),Matrix_I{n+1}(4,5),Matrix_I{n+1}(4,6),Matrix_I{n+1}(5,5),Matrix_I{n+1}(5,6),Matrix_I{n+1}(6,6));
fprintf(fid,'proportional,%f;',0);
n=n+2
end 
system(sprintf('C:\\Users\\User\\Downloads\\documents\\blendyn-master\\MB-Dyn\\mbdyn.exe %s',filename))
n=length(meanx)
close all
Z=importdata('blade10.mov');
if pitch==0
z=Z(n:n:end,4);
figure 
hold on 
time=0:tstep:SimTime+tstep;
plot(time,z)
y=-1.8
plot([0 time(end)],[y y ])
xlabel(' time s')
ylabel(' tip displacement m')
else
    z=Z(n:n:end,4);
figure 
hold on 
time=0:tstep:SimTime+tstep;
plot(time,z)
y=3.4
plot([0 time(end)],[y y ])
xlabel(' time s')
ylabel(' tip displacement m')
end 
    

[Freq,Phase,Resp]=FFTFull(time,z,tstep)
freq=0.54
freq2=0.64
figure 
hold on
plot(Freq,Resp)
plot([freq freq],[0 1])
plot([freq2 freq2],[0 1])
xlabel(' freq Hz')
ylabel(' Magnitude -')