function [Cl_beta]=Cl_beta(CL_wb,Cl_alpha_w,theta,theta_h,Lambda_c2, Lambda_c4,Lambda_c4_h,S_w,b_w,AR_w,AR_h,...
   Z_v,alpha,l_v, M,two_r_one,eta_v,beta,AR_v,b_v,Z_h,x_over_c_v,lambda_v,S_v,...
   Z_w,d,wingloc,Z_w1,S_h,S_o, dihedral,dihedral_h,Lambda_c2_h,l_b,b_h,CL_hb,lambda_w,lambda_h)
% *********************************************
% Cl_beta - Variation of rolling moment coefficient
%        with sideslip angle.
% *********************************************     
%
% Version 8.0 1/23/03 
% This version corrects Cm_alpha, rudder control power, units, and comments.
% AAE565 Spring 2002
%%
% Assumptions:
% Cl_beta is considered to be the sum of the wing-body and horizontal 
%                        tail-body, and vertical tail-body contribution,
% the effect of the fuselage being usually small.
%
% Required Inputs:
% AR_w             Aspect ratio of the wing
% Xw               (positive reaward) distance from the arbitrary reference point to wing aerodynamic center [ft]
% eta_h            Ratio of dynamic pressure at the tail to that of the free stream
% c_w              Mean aerodynamic chord of the wing [ft]
% c_h              Mean aerodynamic chord of the horizontal tail [ft]
% V_h              Horizontal Tail Volume Coefficient
% Xh               Distance from arbitrary refernce point to the horizontal tail aerodynamic center [ft]
% Lambda_c4        Quater Chord Sweep angle [rad]
% B                Compressibility correction factor
% beta             Compressibility correction factor 
% M                Mach number
% Lambda           Mid chord sweep angle
% l_H              Defined in Fig 3.7 of Ref(2) Geometric parameter for horizontal tail location [ft]
% h_H              Defined in Fig 3.7 of Ref(2) Geometric parameter for horizontal tail location [ft]
% Zw               vertical distance from wing root quaterchord point to the fuselage centerline, positve downward
% S_o      Average fusealge crossectional area
%
% A&AE 421 Fall 2001
% Jaret Matthews
% 
% Equations/Figures can be found in :
% Ref. (2) Roskam, Jan. "Methods for Estimating Satbility and
%         Control Derivatives of Conventional Subsonic Airplanes"
%
gama= dihedral;
gama_h=dihedral_h;
Aover_cos_lamc2=AR_w/(cos(Lambda_c2));
Aover_cos_lamc2_h=AR_h/(cos(Lambda_c2_h));
M_cos_lamc2=M*cos(Lambda_c2);
M_cos_lamc2_h=M*cos(Lambda_c2_h);
lf_b=l_b/b_w ;%Length of Fuselage over span, See Ref (2), Pg. 7.15, Figure 7.13
lf_b_h=l_b/b_h ;%Length of Fuselage over span, See Ref (2), Pg. 7.15, Figure 7.13
d=sqrt(S_o/.7854);%Ref. (2) Pg. 7.4 Eq. 7.11
dClB_Gam=-.0005*sqrt(AR_w)*((d/b_w)^2); %Ref. (2) Pg. 7.4 Eq. 7.10 dClB_Gam is in per deg.^2 - Murthy
dClB_Gam_h=-.0005*sqrt(AR_h)*((d/b_h)^2); %Ref. (2) Pg. 7.4 Eq. 7.10 dClB_Gam_h is in per deg.^2 - Murthy
dClB_Zw=(-1.2*(sqrt(AR_w))/57.3)*(Z_w/b_w)*(2*d/b_w);%Ref. (2) Pg. 7.4 Eq. 7.12 - it is in per deg. - Murthy
dClB_Zh=(-1.2*(sqrt(AR_h))/57.3)*(Z_h/b_h)*(2*d/b_h);%Ref. (2) Pg. 7.4 Eq. 7.12 - it is in per deg. - Murthy
% Question Should b_h above be b_w?  Andrisani 2/15/02

%-------------FIGURE 7.11------------------
lam_c_2_fig=[-20 -10 0 10 20 30 40 50 60]*pi/180;
%changed to radians, as input is going to be in radians - Murthy
AR_fig=[1 2 4 6 8];
gam_fig=[1 .5 0];
fig_7_11(:,:,1)=[0.00177	0.00177	0.00177	0.000981	0.000981;...
0.000809	0.000809	0.000809	0.000416	0.000416;...
0	0	0	0	0;...
-0.000949	-0.000949	-0.000949	-0.000714	-0.000478;...
-0.00191	-0.00191	-0.00159	-0.00144	-0.00104;...
-0.00302	-0.00286	-0.00255	-0.00216	-0.00161;...
-0.00429	-0.00405	-0.00366	-0.00311	-0.00233;...
-0.00627	-0.00588	-0.00517	-0.00431	-0.00297;...
-0.00989	-0.00879	-0.00746	-0.00605	-0.00385];

fig_7_11(:,:,2)=[0.00176	0.00176	0.00176	0.00121	0.00121;...
0.00088	0.00088	0.00088	0.000562	0.000562;...
0	0	0	0	0;...
-0.000885	-0.000885	-0.000885	-0.000487	-0.000487;...
-0.00177	-0.00177	-0.00177	-0.00129	-0.00113;...
-0.00273	-0.00273	-0.00241	-0.00209	-0.00185;...
-0.00401	-0.00377	-0.00353	-0.00305	-0.00258;...
-0.00576	-0.0053	-0.00505	-0.00433	-0.00338;...
-0.00887	-0.008	-0.00736	-0.00617	-0.00442];

fig_7_11(:,:,3)=[0.00164	0.00164	0.00164	0.00133	0.00133;...
0.000697	0.000697	0.000697	0.000618	0.000618;...
0	0	0	0	0;...
-0.000882	-0.000882	-0.000882	-0.000645	-0.000645;...
-0.00159	-0.00159	-0.00159	-0.00128	-0.00128;...
-0.00246	-0.00246	-0.00246	-0.00222	-0.00191;...
-0.00341	-0.00341	-0.00333	-0.00301	-0.0027;...
-0.00491	-0.00491	-0.00475	-0.00412	-0.00341;...
-0.00736	-0.00736	-0.00704	-0.00586	-0.00428];

ClB_CL2=interlim3(AR_fig,lam_c_2_fig,gam_fig,fig_7_11,AR_w,Lambda_c2,lambda_w);
ClB_CL2_h=interlim3(AR_fig,lam_c_2_fig,gam_fig,fig_7_11,AR_h,Lambda_c2_h,lambda_h);
%-------------FIGURE 7.12----------------------
A_cos_lamc2_fig=[10	8	6	5	4	3	2];
M_cos_lamc2_fig=[0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9];
fig_7_12=[1	1	1	1	1	1	1;...
1.03	1.02	1.01	1.01	1.01	1.01	1;...
1.06	1.03	1.02	1.01	1.01	1.01	1;...
1.1	1.04	1.04	1.02	1.02	1.01	1;...
1.16	1.07	1.08	1.04	1.04	1.02	1;...
1.23	1.12	1.12	1.08	1.05	1.02	1;...
1.36	1.27	1.18	1.12	1.07	1.02	1;...
1.55	1.41	1.25	1.17	1.1	1.04	1;...
1.84	1.59	1.33	1.2	1.11	1.04	0.992];

KML=interlim2(A_cos_lamc2_fig, M_cos_lamc2_fig, fig_7_12, Aover_cos_lamc2, M_cos_lamc2);
KML_h=interlim2(A_cos_lamc2_fig, M_cos_lamc2_fig, fig_7_12, Aover_cos_lamc2_h, M_cos_lamc2_h);
%-------------FIGURE 7.13-------------------
lam_cos_lamc2_fig=[8 7 6 5.5 5 4.5 4];
lf_b_fig= [0 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6];

fig_7_13=[1	1	1	1	1	1	1;...
0.987	1	1	1	1	1	1;...
0.96	0.979	0.999	1	1	1	1;...
0.92	0.945	0.972	0.993	1	1	1;...
0.87	0.899	0.934	0.966	0.997	1	1;...
0.809	0.845	0.884	0.924	0.972	1	1;...
0.742	0.782	0.828	0.869	0.934	0.978	0.998;...
0.671	0.715	0.765	0.811	0.884	0.946	0.99;...
0.592	0.644	0.696	0.744	0.825	0.908	0.967];

Kf=interlim2(lam_cos_lamc2_fig,lf_b_fig,fig_7_13,Aover_cos_lamc2,lf_b);
Kf_h=interlim2(lam_cos_lamc2_fig,lf_b_fig,fig_7_13,Aover_cos_lamc2_h,lf_b_h);

%-----------FIGURE 7.14-------------------
gam_fig=[1 .5 0];
AR_fig=[1 2 3 4 5 6 7 8];
fig_7_14=[-0.0111	-0.00781	-0.00567;...
-0.00607	-0.00401	-0.00219;...
-0.00369	-0.00227	-0.00101;...
-0.00251	-0.00148	-0.000452;...
-0.00188	-0.00101	-0.000215;...
-0.0014	-0.00061	0.000102;...
-0.00108	-0.000294	0.000339;...
-0.000927	-0.000136	0.000497];

ClB_CL_A=interlim2(gam_fig, AR_fig, fig_7_14, lambda_w, AR_w);
ClB_CL_A_h=interlim2(gam_fig, AR_fig, fig_7_14, lambda_h, AR_h);

%-----------FIGURE 7.15----------------------
Lambda_c2abs=abs(Lambda_c2);
Lambda_c2abs_h=abs(Lambda_c2_h);
lamc2_fig=[0 40 60 ]*pi/180;
%changed to radians, as input is going to be in radians - Murthy

AR_fig=[0 1 2 3 4 5 6 7 8];
gam_fig=[1 .5 0];
fig_7_15(:,:,1)=[0	0	0;...
-5.17E-05	-5.17E-05	-5.17E-05;...
-9.93E-05	-9.93E-05	-8.74E-05;...
-0.000139	-0.000131	-0.000115;...
-0.000169	-0.000157	-0.000133;...
-0.000195	-0.000181	-0.000147;...
-0.000219	-0.000199	-0.000161;...
-0.000236	-0.000213	-0.000169;...
-0.000254	-0.000228	-0.000175];
fig_7_15(:,:,2)=[0	0	0;...
-5.12E-05	-5.12E-05	-5.12E-05;...
-9.54E-05	-9.54E-05	-8.54E-05;...
-0.000134	-0.000126	-0.000112;...
-0.000166	-0.000152	-0.00013;...
-0.000188	-0.000172	-0.000144;...
-0.00021	-0.00019	-0.000156;...
-0.000228	-0.000204	-0.000162;...
-0.000245	-0.000219	-0.000169];
fig_7_15(:,:,3)=[0	0	0;...
-5.75E-05	-5.75E-05	-4.58E-05;...
-9.06E-05	-9.06E-05	-7.31E-05;...
-0.000118	-0.000114	-9.45E-05;...
-0.000139	-0.00013	-0.00011;...
-0.000155	-0.000143	-0.00012;...
-0.00017	-0.000157	-0.000128;...
-0.000182	-0.000167	-0.000133;...
-0.000194	-0.000176	-0.000139];

ClB_Gam =interlim3(lamc2_fig,AR_fig,gam_fig,fig_7_15,Lambda_c2abs,AR_w,lambda_w);
ClB_Gam_h =interlim3(lamc2_fig,AR_fig,gam_fig,fig_7_15,Lambda_c2abs_h,AR_w,lambda_h);

%---------------FIGURE 7.16------------------
A_cos_lamc2_fig=[10 8 6 4 2];
M_cos_lamc2_fig=[0 .1 .2 .3 .4 .5 .6 .7 .8 .9];
fig_7_16=[0.1	0.1	0.1	0.1	0.1;...
0.1	0.1	0.1	0.1	0.1;...
1.01	0.1	0.1	0.1	0.1;...
1.03	1.02	1.01	1.01	0.1;...
1.06	1.05	1.03	1.03	0.1;...
1.09	1.07	1.05	1.04	1.01;...
1.14	1.12	1.1	1.06	1.01;...
1.21	1.14	1.14	1.08	1.02;...
1.32	1.26	1.19	1.12	1.04;...
1.5	1.39	1.28	1.16	1.05];

KMG=interlim2(A_cos_lamc2_fig,M_cos_lamc2_fig,fig_7_16,Aover_cos_lamc2, M_cos_lamc2);
KMG_h=interlim2(A_cos_lamc2_fig,M_cos_lamc2_fig,fig_7_16,Aover_cos_lamc2_h, M_cos_lamc2_h);

%-------------FIGURE 7.17-------------------
AR_fig=[3 4 5 6 7 8 9 10 11];
gam_fig=[1 .6 .4 0];
fig_7_17=[-2.36E-05	-2.36E-05	-2.22E-05	-1.88E-05;...
-3.01E-05	-3.01E-05	-2.87E-05	-2.20E-05;...
-3.34E-05	-3.34E-05	-3.22E-05	-2.37E-05;...
-3.49E-05	-3.49E-05	-3.35E-05	-2.33E-05;...
-3.65E-05	-3.65E-05	-3.39E-05	-2.30E-05;...
-3.66E-05	-3.66E-05	-3.40E-05	-2.40E-05;...
-3.72E-05	-3.72E-05	-3.48E-05	-2.57E-05;...
-3.99E-05	-3.99E-05	-3.67E-05	-2.84E-05;...
-4.00E-05	-4.00E-05	-4.14E-05	-3.25E-05];

dClB_thetatanL4=interlim2(gam_fig,AR_fig,fig_7_17, lambda_w, AR_w);
dClB_thetatanL4_h=interlim2(gam_fig,AR_fig,fig_7_17, lambda_h, AR_h);
%------
% Remember, the values of ClB_CL2, ClB_CL2_h, ClB_CL_A, ClB_CL_A_h are in per deg. and not in per rad. from graph - Murthy
% Remember, the values of ClB_Gam, ClB_Gam_h, dClB_thetatanL4, dClB_thetatanL4_h are in per deg.^2 and not in per rad.^2 from graph - Murthy
% However, as can be seen in equations of the reference, the 57.3 multiplication factor takes care of it
% and finally the results ClB_wb, ClB_hb, ClB_h are in per rad. Only thing we have to make sure is that 
% gama and theta that multiply the ClB_Gam and dClB_thetatanL4 must be in degrees since 
% ClB_Gam and dClB_thetatanL4 are in per deg.^2
ClB_wb=57.3*(CL_wb*(ClB_CL2*KML*Kf+ClB_CL_A)+gama*(180/pi)*(ClB_Gam*KMG+dClB_Gam)+dClB_Zw+theta*(180/pi)*tan(Lambda_c4)*dClB_thetatanL4); %eqn 7.9
ClB_hb=57.3*(CL_hb*(ClB_CL2_h*KML_h*Kf_h+ClB_CL_A_h)+gama_h*(180/pi)*(ClB_Gam_h*KMG_h+dClB_Gam_h)+dClB_Zh+theta_h*(180/pi)*tan(Lambda_c4_h)*dClB_thetatanL4_h); %eqn 7.9
ClB_h=ClB_hb*S_h*b_h/(S_w*b_w); %eqn 7.13
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%The following is used to find Cy_beta_v
%%Some derived constants

b_v_over_two_r_one = b_v/two_r_one;

S_h_over_S_v = S_h/S_v;

Z_h_over_b_v = -abs(Z_h/b_v);

fig7_1_xaxis = 2*Z_w1/d;

%%Interpolating in Figure 7.6
Z_h_over_b_v_fig8 = ...
[-0.014442	-0.075377	-0.17057	-0.30375	-0.38361	-0.43113	-0.48242	-0.52417	-0.56402	-0.6171	-0.67203	-0.71932	-0.74955	-0.7911	-0.82694	-0.87399	-0.9322	-0.97536];
A_v_hb_over_A_v_b_fig8 = ...
[1.2741	1.2141	1.1242	1.0272	0.97877	0.95656	0.94198	0.93866	0.93912	0.95489	0.98205	1.0243	1.0588	1.1161	1.1734	1.2839	1.4552	1.5921]; 
%%x/c_v=.8, Fig 7.6
%data1=[Z_h_over_b_v_fig8',A_v_hb_over_A_v_b_fig8']
A_v_hb_over_A_v_b8 = interlim1(Z_h_over_b_v_fig8, A_v_hb_over_A_v_b_fig8,Z_h_over_b_v);

Z_h_over_b_v_fig7 = ...
[-0.010859	-0.13075	-0.24678	-0.34946	-0.46538	-0.56029	-0.64558	-0.73449	-0.7874	-0.859	-0.89846	-0.9585	-0.97914];
A_v_hb_over_A_v_b_fig7 = ...
[1.2145	1.117	1.0386	0.97908	0.93121	0.92121	0.94924	1.027	1.0889	1.2313	1.3388	1.5306	1.5959]; 
   %%x/c_v=.7, Fig 7.6
A_v_hb_over_A_v_b7 = interlim1(Z_h_over_b_v_fig7, A_v_hb_over_A_v_b_fig7,Z_h_over_b_v);

Z_h_over_b_v_fig6 = ...
[-0.0092336	-0.10051	-0.23366	-0.3439	-0.41987	-0.51294	-0.60968	-0.68163	-0.76667	-0.85531	-0.9117	-0.97913];
A_v_hb_over_A_v_b_fig6 = ...
[1.1381	1.0822	0.99259	0.94081	0.92288	0.89757	0.90669	0.94979	1.0466	1.2007	1.3504	1.5997]; 
   %%x/c_v=.6, Fig 7.6
A_v_hb_over_A_v_b6 = interlim1(Z_h_over_b_v_fig6, A_v_hb_over_A_v_b_fig6,Z_h_over_b_v);

Z_h_over_b_v_fig5 = ...
[-0.01143	-0.16927	-0.27951	-0.40107	-0.532	-0.6552	-0.76115	-0.84601	-0.93625	-0.98294];
A_v_hb_over_A_v_b_fig5 = ...
[1.0541	0.95724	0.90546	0.87295	0.87113	0.9112	1.0006	1.1471	1.3852	1.5959]; 
  %%x/c_v=.5, Fig 7.6
A_v_hb_over_A_v_b5 = interlim1(Z_h_over_b_v_fig5, A_v_hb_over_A_v_b_fig5,Z_h_over_b_v);

Matrixa = [.5, .6, .7, .8];
Matrixb = [A_v_hb_over_A_v_b5, A_v_hb_over_A_v_b6, A_v_hb_over_A_v_b7, A_v_hb_over_A_v_b8];
A_v_hb_over_A_v = interlim1(Matrixa, Matrixb, x_over_c_v);

%%Done interp fig 7.6

%%Interpolating in Figure 7.5
%%The referenced figure assumes all values of lambda_v less than .6 are estimated at .6


if lambda_v < .6;
   lambda_v = .6;
   fprintf('lambda_v is out of limits. reset to 0.6');
end

b_v_over_two_r_one_fig1 = ...
[0.0089225	0.13036	0.2147	0.35595	0.60146	0.79048	1.0743	1.3015	1.5668	1.8701	2.0883	2.2782	2.5345	2.753	3.0286	3.5224	4.0159	4.5378	5.0216	5.5054	5.9986	6.4918	6.9377];
A_v_b_over_A_v_fig1 = ...
[0.10677	0.4613	0.6557	0.85004	1.0595	1.189	1.3374	1.421	1.4741	1.5081	1.5041	1.4696	1.4235	1.3508	1.2666	1.1555	1.0901	1.0629	1.0433	1.0351	1.0346	1.0378	1.0221]; 
   %%lambda_v=.1, figure 7.5
b_v_over_two_r_one_1 = interlim1(b_v_over_two_r_one_fig1,A_v_b_over_A_v_fig1,b_v_over_two_r_one);

b_v_over_two_r_one_fig6 = ...
[0.0079378	0.091375	0.1853	0.27945	0.42102	0.55323	0.7139	1.0262	1.3861	1.6513	1.8789	2.0402	2.2395	2.5245	2.7525	3.0281	3.5316	4.0252	4.5281	5.012	5.5148	6.0081	6.4824	6.9567];
A_v_b_over_A_v_fig6 = ...
[0.29364	0.65965	0.83497	0.96834	1.1017	1.2121	1.3225	1.4633	1.5659	1.6228	1.634	1.63	1.5993	1.5265	1.4462	1.3467	1.2051	1.1359	1.101	1.0738	1.0542	1.0307	1.0264	1.0259]; 
   %%lambda_v<=.6
b_v_over_two_r_one_6 = interlim1(b_v_over_two_r_one_fig6,A_v_b_over_A_v_fig6,b_v_over_two_r_one);

Matrixa = [.6, 1];
Matrixb = [b_v_over_two_r_one_6,b_v_over_two_r_one_1];
A_v_b_over_A_v = interlim1(Matrixa, Matrixb, lambda_v);
%%Done interlim fig 7.5


%%Interpolating in Figure 7.7
S_h_over_S_v_fig = ...
[0.00			0.059518	0.13083	0.19024	0.37989	0.51796	0.66361	0.85212	1.0836	1.3735	1.5652	1.7686	1.9995];
K_H_fig = ...
[0.0039038	0.09466	0.19724	0.28009	0.50117	0.63941	0.75409	0.85712	0.95639	1.0441	1.0763	1.1086	1.1449]; 
  %%Fig 7.7
K_H = interlim1(S_h_over_S_v_fig, K_H_fig, S_h_over_S_v);
%%Done interlim Figure 7.7

A_v_eff = A_v_b_over_A_v*AR_v*(1+K_H*(A_v_hb_over_A_v-1));  %%eqn 7.6

term1 = .724 + 3.06*((S_v/S_w)/(1+cos(Lambda_c4)))+.4*Z_w/d+.009*AR_v; %%eqn 7.5

%%Interpolating in Figure 7.3
b_v_over_two_r_one_fig = [0.04978	0.48968	0.9977	1.506	1.985	2.491	2.9873	3.484	4.0019	4.5002	4.9986	5.5166	5.9855];
k_fig = [0.76234	0.75567	0.76091	0.75441	0.74784	0.83516	0.91853	0.98628	0.98763	0.99283	0.99022	0.98767	0.9928]; 
  %%Fig 7.3
k = interlim1(b_v_over_two_r_one_fig,k_fig,b_v_over_two_r_one);
%%Done interp fig 7.3

CL_a_v = (2*pi*A_v_eff)/(2+sqrt((A_v_eff^2*beta^2/((Cl_alpha_w/(2*pi))^2))*(1+(tan(Lambda_c4))^2/beta^2) + 4));
%%eqn 3.8 with A_v_eff subbed in for A (see pg 7.2)
Cy_beta_v = -k*CL_a_v*term1*S_v/S_w;    %%eqn 7.4
%----------------END Cy_beta_v------------------------

ClB_v=Cy_beta_v*(Z_v*cos(alpha)-l_v*sin(alpha))/b_w;	%eqn 7.14
Cl_beta=ClB_wb+ClB_h+ClB_v; %eqn 7.8
