2016 수학 모델링 인터내셔널 A 프로그램 (오리지널) 저자 cclplus

clear all;
close all;
clc
format long
syms h S Fw Ff Ff1 a b c d l L F depth n pl m x1 y1 y t distance n a1 b1;
F=[];
theta=[]; 
v=24; %  
l=105*10^(-3); %         
L=22.05; %      
num=0; %                    
num1=round(L/l);
num2=0;
lin=0/180*pi; %             
lin1=90/180*pi;
lin2=0;
m2=1200; %     
pg=7.7*10^3; %      （  ：kg/m^3)
depth=20; %  
pl=7; %         
vh=0; %    
g=9.8; %                ，   m/s^2
p=1.025*10^3; %    
M=1000; %    
m=10; %    
m1=100; %        
y=0;
d=1;
j1=0;
j2=0;
while(abs(y-d)>0.005)%           ，
if (y>d)&&(numd)&&(num==round(L/l))
lin1=lin;
lin=(lin1+lin2)/2;
end
%       
Ff1=p*g*pi*(0.3/2)^2;
%       
Ff2=p*g*pi*(0.05/2)^2;
%       
Ffg=p*g*m2/pg;
%          
Fhg=374*pi*((m2/pg/3/4)^(1/3))^2*vh^2;
%            
S=2*(2-h);
%        
Fw=0.625*S*v^2;
%        
V=pi*(2/2)^2*h;
%        
Ff=p*g*V;
%            
Fb=374*2*h*vh^2;
%            
Fs1=374*0.3*vh^2;
%              
Fs=374*0.05*vh^2;
%         
if num==round(L/l)
h=(m2*g+M*g+4*m*g+m1*g-Ff1-4*Ff2-Ffg+pl*L*g+(Fhg+4*Fs+Fs1)*tan(lin))/(p*g*pi-(1.25*v^2+374*vh^2)*tan(lin));
else 
h=(m2*g+M*g+4*m*g+m1*g-Ff1-4*Ff2-Ffg+num*pl*l*g)/(p*g*pi);
end
a=Fw+Fb;
b=-M*g+Ff+(Fw+Fb)*tan(lin);
if j1==0
a=eval(a);
b=eval(b);
else
end
F(1)=sqrt(a^2+b^2);
theta(1)=atan(b/a);
n=0;
for i=1:4
%          
Fh(i)=374*0.05*sin(theta(i));
n=n+Fh(i);
a=Fw+Fb+n;
if j1==0
a=eval(a);
else
end
b=F(i)*sin(theta(i))+p*g*pi*(50*10^(-3)/2)^2-m*g;
F(i+1)=sqrt(a^2+b^2);
theta(i+1)=atan(b/a);
end
c=0;
for i=1:5
c=c+sin(theta(i));
end
d=depth-c-h;
y1=lin;
distance=0;
if num==round(L/l)
y=l*sin(y1);
x1=Fw/sqrt(1-(sin(y1))^2);
for i=1:num-1
m=(x1*sin(y1)+i*pl*l*g)/sqrt((x1*sin(y1)+i*pl*l*g)^2+Fw^2);
m=m*l;
y=y+m;
n=Fw/sqrt((x1*y1+i*pl*l*g)^2+Fw^2)*l;
if j1==0
n=eval(n);
else
end
distance=distance+n;
if j1==0
y=eval(y);
else
end
end
else
y=y1*l;
distance=(round(L/l)-num)*l;
for i=1:num
x1=Fw/sqrt(1-(sin(y1))^2);
m=(x1*y1+i*pl*l*g)/sqrt((x1*y1+i*pl*l*g)^2+Fw^2)*l;
y=y+m; 
n=Fw/sqrt((x1*y1+i*pl*l*g)^2+Fw^2)*l;
if j1==0
y=eval(y);
n=eval(n);
else
end
distance=distance+n;
end
end
m=0;
j1=1;
j2=j2+1;
end
%       
Ff1=p*g*pi*(0.3/2)^2;
%       
Ff2=p*g*pi*(0.05/2)^2;
%       
Ffg=p*g*m2/pg;
%          
Fhg=374*pi*((m2/pg/3/4)^(1/3))^2*vh^2;
%            
S=2*(2-h);
%        
Fw=0.625*S*v^2;
%        
V=pi*(2/2)^2*h;
%        
Ff=p*g*V;
%            
Fb=374*2*h*vh^2;
%            
Fs1=374*0.3*vh^2;
%              
Fs=374*0.05*vh^2;
%         
if num==round(L/l)
h=(m2*g+M*g+4*m*g+m1*g-Ff1-4*Ff2-Ffg+pl*L*g+(Fhg+4*Fs+Fs1)*tan(lin))/(p*g*pi-(1.25*v^2+374*vh^2)*tan(lin));
else 
h=(m2*g+M*g+4*m*g+m1*g-Ff1-4*Ff2-Ffg+num*pl*l*g)/(p*g*pi);
end
a=Fw+Fb;
b=-M*g+Ff+(Fw+Fb)*tan(lin);
F(1)=sqrt(a^2+b^2);
theta(1)=atan(b/a);
n=0;
for i=1:4
%          
Fh(i)=374*0.05*sin(theta(i));
n=n+Fh(i);
a=Fw+Fb+n;
b=F(i)*sin(theta(i))+p*g*pi*(50*10^(-3)/2)^2-m*g;
F(i+1)=sqrt(a^2+b^2);
theta(i+1)=atan(b/a);
end
disp('            （   ）')
th=90-theta*180/pi
m=85*pi/180;
if theta(5)>m
disp('         ，    ')
else 
disp('        ')
end
c=0;
for i=1:5
c=c+sin(theta(i));
end
d=depth-c-h;
y1=lin;
distance=0;
if num==round(L/l)
y=l*sin(y1);
x1=Fw/sqrt(1-(sin(y1))^2);
for i=1:num-1
m=(x1*sin(y1)+i*pl*l*g)/sqrt((x1*sin(y1)+i*pl*l*g)^2+Fw^2);
m=m*l;
y=y+m;
n=Fw/sqrt((x1*y1+i*pl*l*g)^2+Fw^2)*l;
distance=distance+n;
plot(distance,y,'o')
hold on
end
else
y=y1*l;
for i=1:round(L/l)-num
distance=i*l;
y=0;
plot(distance,y,'o')
hold on
grid on
end
for i=1:num
x1=Fw/sqrt(1-(sin(y1))^2);
m=(x1*y1+i*pl*l*g)/sqrt((x1*y1+i*pl*l*g)^2+Fw^2)*l;
y=y+m; 
n=Fw/sqrt((x1*y1+i*pl*l*g)^2+Fw^2)*l;
if j1==0
y=eval(y);
n=eval(n);
else
end
distance=distance+n;
plot(distance,y,'o')
hold on
end
end
m=0;
for i=1:5
m=m+cos(theta(i));
end
%       
disp('         ')
ans=distance+m