I am unable to understand how to execute this code.

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

Answer 1

[James Tursa]

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.)