I am unable to understand how to execute this code.

3 次查看(过去 30 天)
clear all
n = input('enter the value of no. of day');
% t = input('enter the value of time');
H = input ('enter the altitude');
CF = input ('enter the fraction of cloud present');
V_wind = input ('enter the value of velocity of wind');
% clear t
%%Constants taken
u = [7.447323 2.071808 9.010095 7.981491 194];
a = u(1); % 7.447323;
b = u(2); % 2.071808;
c = u(3); % 9.010095;
d = u(4); % 7.981491;
l = u(5); % 194;
%%Profile
x = (0:0.1:u(5));
f = 1/8.*(sqrt(u(1).*(u(5)-x).*(u(2).*x-u(5)*sqrt(u(3))+sqrt(u(3)*u(5).^2-u(4)*u(5).*x))));
D = 2*max(f);
R = D/2;
clear x
%%Area of Surface
syms x
y = 1/8*(sqrt(u(1)*(u(5)-x).*(u(2).*x-u(5)*sqrt(u(3))+sqrt(u(3)*u(5).^2-u(4)*u(5).*x)))); % Given curve for shape
y1 = diff(y);
y2 = sqrt(1+y1.^2);
y3 = y*y2;
A = 2*pi*int(y3,0,u(5));
vpa(A,5) % Area of Envelope
A = ans
m = 120*10.^-3*A*0.12*10.^-3
clear x y
%%constants
alpha = 0.15;
albedo_ground = 0.35;
albedo_cloud = 0.5;
P_sea = 101325;
e = 0.0167;
%%solar model
syms t
delta = 23.45.*sin((pi/180).*(360/365).*(284+n));
% [T_oa, P_oa, rho_oa] = atmos(H);
P_oa = 5446;
rho_oa = 0.088;
T_oa = 216.65;
w = 15.*(12-t);
phi = 18.975;
theta = asin(sin((pi/180).*delta).*sin((pi/180).*phi)+cos((pi/180).*delta).*cos((pi/180).*phi).*cos((pi/180).*w));
% psi = asin((sin(pi/180).*w).*cos((pi/180).*delta)/cos(theta));
MA = 2.*pi.*n/365;
M = (P_oa/P_sea).*(sqrt(1229+(614.*sin(theta)).^2)-614.*sin(theta)); %call fun
Tau_atm = 0.5.*(exp(-0.65.*M)+exp(-0.95.*M));
zita = MA+2.*e.*sin(M)+1.25.*e.^2.*sin(2.*M);
ID = (1367.*Tau_atm.^M).*((1+e.*cos(zita))./(1-e.^2)).^2;
I_d = (1/2).*ID.*sin(theta).*(1.-Tau_atm.^M)./(1.-1.4.*log(Tau_atm));
Q_D = alpha.*ID.*sin(theta) % direct radiation
% diffuse soalr radiation model
Q_d = alpha.*I_d
%Reflected solar radiation model
albedo = albedo_ground.*(1-CF).^2+albedo_cloud.*CF;
Ir = albedo.*(ID.*sin(theta)+I_d);
omega = pi.*abs(12-t)/12;
Q_r = alpha.*Ir.*cos(omega)
%Infrared Radiation from earth
T_g = 280; %temperature of ground in k
T_cl = 262.15; % assuming cloud at altitude 4km in k
R = 6371000; %redius of earth
epsilon_e = 0.95;
alpha_ir = 0.85;
sigma = 5.67.*10.^-8;
tau_atm_IR = 1.716-0.5.*(exp(-0.65.*P_oa/P_sea) + exp(-0.95.*P_oa/P_sea));
T_e = T_g.*(1-CF)+T_cl.*CF;
eta = asin(R/(R+H));
VF = (1-cos(eta))/2;
Q_earth = sigma.*epsilon_e.*alpha_ir.*VF.*tau_atm_IR.*T_e.^4
%Infrared Radiation from sky
Pw = 0.0042;
epsilon_sky = 1;
T_sky = T_oa.*(0.51+0.208.*sqrt(Pw)).^0.25;
Q_sky = sigma.*epsilon_sky.*alpha_ir.*(1-VF).*T_sky.^4
%Infrared Radiation from airship
syms T_f
epsilon_f = 0.85;
Q_a = 2.*sigma.*epsilon_f.*T_f.^4
r = 0.05;
Q_aIR = alpha_ir.*sigma.*epsilon_f.*T_f.^4/(1-r)
% Convection heat transfer model
% external convection
g = 9.81;
K_air = 0.0241.*(T_oa/273.15).^0.9;
beta_air = 1/T_oa;
Pr_air = 0.804-3.25.*10.^(-4).*T_oa;
nu_air= 1.46.*10.^-6.*T_oa.^1.5/(rho_oa.*(T_oa+110.4));
Ra_air = g.*beta_air.*(T_f-T_oa).*D.^3/nu_air.^2.*Pr_air;
Re_air = V_wind.*l/nu_air;
h_free = (K_air/D).*(0.6+0.387.*((Ra_air)/(1+(0.559/Pr_air).^(9/16)).^(16/9)).^(1/6)).^2;
h_forced = (K_air/l).*(Re_air.*Pr_air.*(0.2275/(log(Re_air)).^2.584-850/Re_air));
h_oc = (h_free.^3 + h_forced.^3).^(1/3);
Q_oc = h_oc.*(T_f-T_oa)
% %internal convection
% syms T_he
% rho_he = 0.1786
% k_he = 0.144.*(T_he/273.15).^0.7
% Pr_he = 0.729-1.6.*10.^(-4).*T_he
% beta_he = 1/T_he
% nu_he = 4.32.*10.^-6.*(T_he.^0.674)/rho_he
% Ra_he = g.*beta_he.*(T_f-T_he).*D.^3/nu_he.^2.*Pr_he
% if Ra_he<1.5.*10.^8
% Nu_ic = 2.5.*(2+0.6.*Ra_he.^(1/4))
% else Nu_ic = 0.325.*Ra_he.^(1/3)
% end
%
% h_ic = Nu_ic.*k_he/D
% Q_ic=h_ic.*(T_f-T_he)
%%solution
C_p = 1.42;
h = 1;
x = 0:h:23
y = zeros(1,length(x))
y(1) = 216.65
F_xy = (Q_D+Q_d+Q_r+Q_earth+Q_sky-Q_aIR-Q_a-Q_oc)/(m.*C_p)
for i=1:(length(t)) % calculation loop
k_1 = F_xy(x(i),y(i))
k_2 = F_xy(x(i)+0.5.*h,y(i)+0.5.*h.*k_1)
k_3 = F_xy((x(i)+0.5.*h),(y(i)+0.5.*h.*k_2))
k_4 = F_xy((x(i)+h),(y(i)+k_3.*h))
y(i+1) = y(i) + (1/6).*(k_1+2.*k_2+2.*k_3+k_4).*h % main equation
end

采纳的回答

James Tursa
James Tursa 2015-6-8
编辑:James Tursa 2015-6-8
To answer your question literally, put the above code in a file with the extension .m, and then type the file name without the extension. E.g., if you put this in a file called mycode.m, then type at the command line: mycode
If you want more help than this, then please provide a more detailed question (what problems you are having when running the code, any errors you are getting, etc.)

更多回答(0 个)

类别

Help CenterFile Exchange 中查找有关 Parallel Computing 的更多信息

标签

尚未输入任何标签。

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by