From 83a23144093472736bfb417a31a764e914234fc4 Mon Sep 17 00:00:00 2001 From: Charles N Wyble Date: Tue, 10 Dec 2024 11:28:49 -0600 Subject: [PATCH] capture to git --- Balloon_Shape/create_balloon.m | 221 +++++++++++ Balloon_Shape/natural_balloon.slx | Bin 0 -> 13764 bytes Common_Functions/atmo_p.m | 61 +++ Common_Functions/color_line.m | 43 +++ Common_Functions/coordinates_to_dxdy.m | 26 ++ Common_Functions/dxdy_to_coordinates.m | 26 ++ Common_Functions/plotxx.m | 87 +++++ Common_Functions/savex.m | 143 +++++++ Common_Functions/sphere_CD.m | 8 + Common_Functions/stdatmo.m | 405 ++++++++++++++++++++ Common_Functions/stdatmo_H2.m | 22 ++ Common_Functions/transpose_arrayOfStructs.m | 16 + Common_Functions/volRevolve.m | 157 ++++++++ Common_Functions/wraptopi.m | 28 ++ GUI.fig | Bin 0 -> 5729 bytes GUI.m | 251 ++++++++++++ GUI_inputs/default_GUI_inputs.mat | Bin 0 -> 330 bytes balloon_inputs.mat | Bin 0 -> 455 bytes license.txt | 24 ++ 19 files changed, 1518 insertions(+) create mode 100644 Balloon_Shape/create_balloon.m create mode 100644 Balloon_Shape/natural_balloon.slx create mode 100644 Common_Functions/atmo_p.m create mode 100644 Common_Functions/color_line.m create mode 100644 Common_Functions/coordinates_to_dxdy.m create mode 100644 Common_Functions/dxdy_to_coordinates.m create mode 100644 Common_Functions/plotxx.m create mode 100644 Common_Functions/savex.m create mode 100644 Common_Functions/sphere_CD.m create mode 100644 Common_Functions/stdatmo.m create mode 100644 Common_Functions/stdatmo_H2.m create mode 100644 Common_Functions/transpose_arrayOfStructs.m create mode 100644 Common_Functions/volRevolve.m create mode 100644 Common_Functions/wraptopi.m create mode 100644 GUI.fig create mode 100644 GUI.m create mode 100644 GUI_inputs/default_GUI_inputs.mat create mode 100644 balloon_inputs.mat create mode 100644 license.txt diff --git a/Balloon_Shape/create_balloon.m b/Balloon_Shape/create_balloon.m new file mode 100644 index 0000000..a93e8c6 --- /dev/null +++ b/Balloon_Shape/create_balloon.m @@ -0,0 +1,221 @@ +function balloon_output = create_balloon(balloon_input_filename) +% Authors: +% Brent Justice - justiceb@purdue.edu +% Ian Means +% Purpose: +% Design contour for zero pressure balloon +% With given mission criteria +% References: +% "A Parallel Shooting Method for Determining the Natural +% shape of a Large Scientific Balloon", Frank Baginski, +% William Collier, Tami Williams +% +% "Scientific Ballooning", Nobuyuki Yajima, Naoki Izutsu, +% Takeshi Imamura, Toyoo Abe +% Units: +% ALL calculations performaed in SI units. Only converted +% to english units for printing or plotting + +%% Inputs (Mission Criteria) +load(balloon_input_filename) %load the following input parameters +%rho_PE --> kg/m^3 density of polyethlyene we purchased +%thickness_PE --> m thickness of polyethylene we purchased +%Wpayload --> N payload weight +%alt_apogee --> altitude at apogee (m) +%numGores --> number of gores +%R_gas --> gas constant of lifting gas + +%constants (mission criteria and materials) +g = 9.81; %m/s/s +R_air = 287.058; %specific gas constant air (SI) +%R_H2 = 4124; %specific gas constant hydrogen (SI) +%R_He = 2077; %specific gas constant helium (SI) + +%% float conditions (at apogee) +[rho_air,a_apogee,T,P,nu,ZorH]=stdatmo(alt_apogee); %SI (standard atmosphere) +rho_gas = P/(R_gas*T); %Ideal gas law, assume ambient pressure +a = 0; %distance from zero pressure line to base (0 for apogee) +k = (2*pi)^(-1/3); %constant of geometry +tb = 1; % == b/bd == 1 at apogee +rho = 1; % == w/wd == 1 at apogee +bd = (rho_air - rho_gas) * g; %effective bouyant force per unit volume at apogee (N/m^3) +wd = rho_PE * thickness_PE; %mass per unit area of PE at apogee (kg/m^2) [changes with altitude] +lambda = (Wpayload/bd)^(1/3); %shaping parameter (m) +epsilon = (2*pi)^(1/3) * wd * g / (bd*lambda); %shaping parameter [] +a_dash = a/lambda; + +%print sizing parameters in english units (easier to digest) +%fprintf('bd = %f lbs/ft^3 \n',bd / 157.087464) +%fprintf('Payload = %f lbs \n',Wpayload / 4.44822162) +%fprintf('lambda = %f ft\n',lambda*3.28084) +%fprintf('wd = %f lbs/ft^2 \n',wd *0.204816144) +%fprintf('epsilon = %f \n',epsilon) +%fprintf('Number of gores = %f\n', numGores) + +%% Shooting method simulation for finding gore contour shape +%constant initial conditions +r0 = 0; %radius at base = 0 +z0 = 0; %height at base = 0 +S0 = 0; %area at base = 0 +V0 = 0; %volume at base = 0 + +%final value conditions +r_end_condition = 0; %radius at apex = 0 +theta_end_condition = -pi/2; %angle at apex = -90 (radians) + +options = simset('SrcWorkspace','current'); +%iterate +theta_end_error = 0.1; %set error high to enter loop +theta_0_guess = 60 * 0.0174532925; %angle at base (rads) [GUESS] +m0 = 2*pi*cos(theta_0_guess); %calculation dependant on guess +while abs(theta_end_error) > 0.00174532925 %while error high, keep iterating + r_end_error = 2; %set error high to enter loop + s_dash_end_guess = 2; %arclength of gore (m/lambda) [GUESS] + while abs(r_end_error) > 0.001 %shooting method until error for r is small + sim('natural_balloon',[],options) %run simulink model + s_dash = simout.time; %format sim outputs + z_dash = simout.data(:,1); + r_dash = simout.data(:,2); + S_dash = simout.data(:,4); + V_dash = simout.data(:,5); + theta = simout.data(:,3); + + %Calcuate non-dashed values + r = r_dash*lambda; %m + s = s_dash*lambda; %m + z = z_dash*lambda; %m + S = S_dash(end)*lambda^2; %balloon surface area (m^2) + V = V_dash(end)*lambda^3; %balloon volume (m^3) + + %calculate radius error + r_end_error = (r_end_condition - r(end)); + theta_end_error = theta_end_condition - theta(end); + + %Use interpolation to determine next guess + if exist('shoot1') == 0 %this was our first guess. We must pick a 2nd arbitrary guess + shoot1.A2 = s_dash_end_guess; + shoot1.Z2 = r(end); + + s_dash_end_guess = s_dash_end_guess + 0.0001; %set next guess to our first guess + some small amount + else + shoot1.A1 = shoot1.A2; + shoot1.Z1 = shoot1.Z2; + shoot1.A2 = s_dash_end_guess; + shoot1.Z2 = r(end); + Zsol = r_end_condition; + a = Zsol - shoot1.Z1; + b = shoot1.Z2 - Zsol; + + Anew = (a*shoot1.A2 + b*shoot1.A1)/(a+b); %calculate next guess based upon interpolation of previous 2 results + s_dash_end_guess = Anew; + end + end + clearvars shoot1; + + %calculate theta_0 guess for next iteration using interpolation method + if exist('shoot2') == 0 + shoot2.A2 = theta_0_guess; + shoot2.Z2 = theta(end); + + theta_0_guess = theta_0_guess + 0.001; + else + shoot2.A1 = shoot2.A2; + shoot2.Z1 = shoot2.Z2; + shoot2.A2 = theta_0_guess; + shoot2.Z2 = theta(end); + Zsol = theta_end_condition; + a = Zsol - shoot2.Z1; + b = shoot2.Z2 - Zsol; + + Anew = (a*shoot2.A2 + b*shoot2.A1)/(a+b); + theta_0_guess = Anew; + m0 = 2*pi*cos(theta_0_guess); + end +end + +%% Calculate sea level hydrogen needed +m_gas = rho_gas*V; %kg --> mass of lifting gas needed +[rho_air_SL,a_SL,T_SL,P_SL,nu_SL,ZorH_SL]=stdatmo(0); %SI (standard atmosphere) +rho_gas_SL = P_SL/(R_gas*T_SL); %Ideal gas law, assume ambient pressure +V_gas_SL = m_gas / rho_gas_SL; %m^3 lifting gas Volume at Sea Level + +%% Print outputs +m_PE = S*wd; %kg +m_payload = Wpayload / g; +arclength = s(end); +height = z(end); +fprintf('epsilon = %f \n',epsilon) +fprintf('balloon surface area = %f m^2 \n',S ) +fprintf('balloon volume = %f m^3 \n',V ) +fprintf('balloon volume = %f ft^3 \n',V *35.3147) +fprintf('balloon mass = %f kg \n',m_PE ) +fprintf('payload mass = %f kg \n',Wpayload/g ) +fprintf('total weight = %f N \n',Wpayload+m_PE*g ) +fprintf('total lift = %f N \n',V*bd ) +fprintf('lifting gas Volume at Sea Level = %f ft^3 \n',V_gas_SL *35.3147 ) +fprintf('angle at apex = %f degrees (shoot for -90) \n',theta(end)/ 0.0174532925 ) +fprintf('angle at base = %f degrees (shoot for -90) \n',theta(1)/ 0.0174532925 ) +fprintf('radius at apex = %f m (shoot for 0) \n',r(end) ) +fprintf('arclength of gore = %f m \n',arclength ) +fprintf('\n') + +f1 = figure(1); +subplot(1,2,1) +plot(r,z) +xlabel('radius to center of balloon (m)') +ylabel('distance from balloon base (m)') +grid on +axis equal +title('balloon gore contour at apogee') +hold all + +%% Determine Gore Shape +goreTheta = (2*pi)/numGores; %radians rotation per gore +s_print = 0:0.1:s(end); +for n = 1:1:length(s_print); + r_print = interp1(s,r,s_print(n)); + goreWidth(n) = r_print*goreTheta; %m arc length formula + fprintf('At %.1f cm from bottom, gore is %.2f cm from center\n',s_print(n)*100,goreWidth(n)/2*100) +end + +subplot(1,2,2) +plot(goreWidth/2,s_print,-goreWidth/2,s_print) +xlabel('width (m)') +ylabel('length (m)') +title('Flat Gore Shape (stencil)') +axis equal +grid on +hold all + +%% 3D plot +y = r; +x = z; +ni = length(x); +nj = 2*ni-1; +tt = linspace(0,2*pi,nj); +for i=1:ni +for j=1:nj + xx(i,j) = x(i); + rr(i,j) = y(i); + yy(i,j) = rr(i,j)*cos(tt(j)); + zz(i,j) = rr(i,j)*sin(tt(j)); +end +end + +f2 = figure(2); +surfl(yy,zz,xx) +colormap(gray) +axis('equal') +xlabel('(meters)') +ylabel('(meters)') +zlabel('height (meters)') +title('Balloon Shape at Apogee') +axis equal + +%% Output +w = whos; +for a = 1:length(w) +balloon_output.(w(a).name) = eval(w(a).name); +end + +end diff --git a/Balloon_Shape/natural_balloon.slx b/Balloon_Shape/natural_balloon.slx new file mode 100644 index 0000000000000000000000000000000000000000..3149841dd7c630fa18d0a330947cf94f01a8b3a7 GIT binary patch literal 13764 zcmaL81B@s^_x3rqZQHhO+q`4jwr$(CZQHhW$1{81EWZ7}*>Ag(u5@?$m&#M8I&~`N zQIG}(K>+{&fB-l*$rG6ZRKzj|1OV6p0sz4N_teJJ+0ew$*^u7Y&e2rC(azq~(b>Y( ziO$`|Iz&lZevkp_~gbH+u9XvDNYY)fE_wBEoKqv+Rs} zuO$0cCCR$q!{QkvfkqtC+%|am1*C`m&jki49|Nt1F^kMR9*ICh{=s*-Rnoj(F zl6f)jiEH30p+;oGnU%t{nPLHP2a{X_f}d^-?U8;p7|-|b<0(M(nGU9yl8xmDq38fnXW5O`dGs4#}P1IP~2#jh_Hq27D1au<- zu-3LR(+AB%-D}Ok!B)E^WtLcaC+dsF6!iVJB=%utRctqqQFrJ&LtK}>Bm(?r)0}?g zTwUyXA1urgk$ybl%$T#i+d}_DKJ>l;@s@P;pD=jmHdDivAW<|B9XQZD2>71fxH6&Z^ zKau{n<6tO#nD=@_0x7in-J(box#$oe2$f7v4TYr2;V`Z;5_i#JV%g4xxC=o+wam0NJt@7DtD5x4pK4rj=sQ@-jW zH#Q?%I9D*)yClck%LR!v)h7|W5 z78Z6!=NOT>PjgT2|&PFCUi>AxT%SQWINC zcPX=M6bzo-SYLlnXUT14+7YNO8X~Ttp>g>Q%nKX@JaydO9>UzbCbh>p51)Du11Gi` z9v)uB&4QVcfs_j4F@^e6Cho^ktGjgnqS{MK?8Nr><;&eTM#PBY&|krd{Px-%EiI)( zX65+AqUbQQvg#c>tCgILEV`NtIJ!AkLp$yRj{aGSt@`!+Pz2U}K5!jWa}owzl?GoKML_t8#|V6(@L6Lv_gKfI12&BAzwdD!IUP?ZIT8;`TSXermCCb)c#=jQ4JzlJaSTzY_^pt1DQW_Gf3S}bm-zxYgsUl6vk5^~`vX?`||ipCz57>{b_ z(Ao+U6Vr|-r2Mcf{ZQaQxR{uqKB(#I;#gU&`=xSpbZ#>;8Wt88fdoC9y1%ORWe%Y{qh;^DAE3Q;WhUc@=(rIaFIXf;6L&L&?GC*NZ0Tn}KdLcpvv$wZLsUG8Q z0?=bv&^^O$_U+%uzt^~;YT3L|GO)tY3=#|H@=8jI-5<`1PoZ&*kK8|JJ0uIbA|sny zmWH($1uQ!o*6vW2{DjNRD-2hrP*Nlcxw;l z^4l$P(z3Yu;|rYEoE~Nhw~voW^2^ZjII9VMzO)(IUoM8BJN|6(H4n$f4ZFi*ro_kR zH}rAWz+WiI&W?w}LB5OA<8kqZ@N)mP+`d7qaUKR7#Kgp&s;$Jztp!E8x`;xLFD@n_ zK(V2w1|lLNqH(fPRaNc(a8=FnjSL9%yn zkQ{+gkm-6oIeB<}tyEA{O!0enjZo+S!^_T& zZN6dzp~UJrn?i&>L>u!2urEInm@J0`spsM3qNXMjS;S6Xz$jr!!T0B7i#FiG%gbxE zhpm@;^VKLiKGY+o0G%A!MQE9}UDOifre`?p4*aP$5aoAZ9p_H)7#?0B z+T;q=ZJJ~=x}q~Y%shUmh#>9bfZdfrs!+J=-tonkcfhdQoyb=?C8z?jABAUGYPFxRNIk7^(CWA^(5ff-Ckc$k@jd=sre1VmcEiQXOWT7+wHa? zC;=@AWyR<7Ma^<^bDptrK}t#+J_gpQiD|&`;yD1=tXgA3lb3CD#MExblUCbR?V4Sx z5wl{Dj42j(H2#8Vw=?0885^8eMPp?ooT*@}!{jpt6bVc3+tb;cg1$Zq3K~cG=q=2$ z@Xw}&gSUlR>Ib98EPq|0Q0Ea1SBmuG*kyTKWoFH6KXL32?=y^10&-=gue66p&DYILL=9ORQAKeEChXp70z`0sPfIo=1v=cklvE!uz`86w&GFeajn9qq^L+fn;^n2X zf4hJJAi}U21>(LyYInReHM*U+N^|>4CrS_Pz2;XWq1O|{+er(fD9KLdVTNhbj|75&BcZ0vd)S>VH#iGtHY>`w+N0LiVBj1eJHg? ze#;>*VwAW3)g_0&C(9P*`fE{mOs^K$7>oKTal-bMEz0K3PR-z#XYf9!FlB4spC*ZdS3z9t( z4Gm4Xp0xJXLpMh1Gs@<@ct$#HNZ3+{9>-=$#XT^U6o`3!XBD^LeN|QPNZO_4r`a6W zb{jDk77C^AEwDpzcQj)8lLiD<`_n>s(as9lmMwT9F3zYu@WI9VxO9alI4D?<&v`Qw zz&tVp$ad&1yVr~QiHi&9$Hy2o^*jADDrJU>JhA3G%+_@ds0_ajM^G$adX{i8ikG7CfKJA6NIzyZ)w!y*|-Dj3gV+CE>K2nBSj8=SGTihTzQJvg`}nL4_o%|G5|adjaf zgU~ApRe7?oFh-Us{IKX)kI>M9gBO596jvV;pQXR=7gi~SUXc(k;mM7TDEl)&G}`oX z5|#*_A5y+41eAJImUsCDJGv%Dj!bC_v|4OZr3QV-g*+CeT6JjQ*Xy%mRvgwA%1bWp z88bhjAt?z7y6SuCxJ98wk(Zb8$V?H-_J53R_38`Uz!s#3acupzw?TQz>G#p6W=1IR zVi}^GmVyBHZ&~!SVLS1^StEVNW|cgaq@*~~54bLiH(C*0 zQ0?7O@tjIZkOWo}=+c}=@0$U8IV#&VtP@9tK&1=EfuNd{UY&4psaUXh?;P3|dIxo? z^Xu=oc5`XWf5#UF78b&8?Qs~w0s>IfWP8`N;JQ8i5ny3qH&oM^8=Jy>?zuyKD59+5 zQ`l`^tJ|1H4esmU*z#Dic%^SmS!oFgd%7=c4$Za```la~hg|7XfpHmGJ0QJ-)+VK- zbjF54%+1Zwpbn+~t)VBgGK~PB1+|i9m-F}RGDXYrw^^*^02;}!p3Km-fyMG8m#uXU zA=i!QdS<*5ui2N0bD-co$voEK#gr+BN3BOZ>PliaS} z)+;))Lq5|X|nh@a2GU;JlE zkU=i9w>4suWt@l2BubOgNf)qsof-pTQ1p>VCki6SQ}svjkl}&lJ#N8A4KoAz2T_-c z`fp5bbZPWk9KU_sHuB^;@#l$1D#k1yWDZ ztRk7g(|XFHX0Ss#;-G=6r(TLjkM%_c450(D$?G>5VF#=D5_ySx$%ukDJUZODrdi<` zN2LXyeCV+Qt^LCI-6-U@h`=sB(AYF^~t?^<~ zY5jO_XPrN2B&4an7fQsz$lkHy?iJI8Qzdz{f zzD#Y>6a40dVnf%@Em9o-E!Pd{*P}B^V)6H>*YE8hA0Dg4&wZX72M#{;j_WQb=qIcI zcn{_vYIk@4d|%tCviQB+SMm3W1}e>&)kWA%RO63Mmfq9TsvGiS#3(~9dj04Lf1Ur7 z%!QSp)kPxG@AFAu80k*Z?+2xXFJ_74^AgqjNI0-uaPwFutu3?2v**+Bj&&FL_wMe% zz`nVi{+cqLmI*H>{$np$!5Bm+<*z;1rZ4EH(&$XlvC-=>71UJjX;H+)m=NJ}b>5Hy zP+$>Plt9!~trjFF)DXo13uZt41+UWvw0ko+1bAxzH1btcd%rr7Kko6*oX8-SH=BQS zKsW>2Yg|1$0>ht!Ku~8HBZ$>Mw(X5b2%=kWTs&ql$#qD-$u%=nAW800K9d5t{z~lT zjIEq9s4_(AO3WsBCUQIjpsbkPig>}Z5;&lgN=9P>SKzVpFyPgsKdeuRL@vp>{K#>y zkVKjBW^&a&uDW8s$BU{d%&fRP^CNNIxfg1C-iS!XsR(UB;0^Z}J4CQg!2wX33uuZc z*`1yf07nVM4Uj@$M1i;+NSGi(ME(YR-f)IOrMBAu{6FQ*AYHp6&_w8s)p6sC9r^5Rib3l={uNKVHdBw~8S~*m3il?xiDl9y?&kWq@TkgRZ$$GcH7v$|VShZ;qqKw& z)K9_OUTqOXGKSGO5{f~!=@gVP9dCVR_r;?-=t~X+o>?Q$v=1fGXDOo!K$}@*6{t?9T=0o%0DrOR1z{h7l?;-YT@jq`rkY-Gxo<n4S@V!cchiK@zNzK_Gh!oJW+I85#{MA4Etss@A9mcJDdutqslq z(e3Tw529dQX@K3$Ta|Q~h0&+&33clYeIkMKtb%)d28A8CbNT6_xO{87h2a46N^gG8qF~`7qv{@h=&N_=Y=7V zBlzgJL*bBEi(!pIWGWZnsYjiM`XdYNQPH%0dlc*@3sT&Zt-)sfDL!sbbQ_8jm46{9 zgu^sT7$S!NYB(fEyrH(FTG`^`oIvkXSDe6+otDog1FQWo~9BW#4jU^)3HHRx(hOI^9m7x2+$oMfGyDu0X8~a2ALUuENp}?Vwogpwl>!tDl_w26kHVqZ zTLLW1unnu0u|SRzP-N->3W4Dzh!D?9RgWB90>FW>81#TsL-Ls59Kn{w_c(03B#?RC zJbJB(8@Ic9D;~TpLjX-T0?W^6v0&-r)pm;Koe3l-2YF`Djak_uh-oXy<`gRparMb2 z50GfH&l~-{7lhQ(Dg1#k`p|m!O9B2K#e_s&)y#W9L=9*~+}Yz`3F`wPgw*}TNzIiL zk0r>r^u8{g+px0uw9m%YEUeY)_vE1UDURD44$aIDFY{ASA(JnnAh5-VJ5HwK_0YhQM402hB2Fp_7s4C z_NGb&3PMLTLc09}(`$tkIGpH>pEgS~z+uSd&#O+70lu(m5)GJ|mRf2t57#*=2wb1< zip19mbGSsZ_gyMFUM%?M!ts5n-%53$0CT2Jcl_d8dCjXG0JR*gAPX($;alw_V|5^x z-%hY*IM8>sR+DPd)6SeDXg~p*Sd@$eQeE*Dtu!0dd2Xv0!wCP!ixPNe&BRNpX^9u` z7C>D|jEkl!WnF`rY|CUOh6%$J|#GM{2=euHYbZZRux^; z;Z>knw@LV8WI*rC$q>lMFa~D+&sKpzj(x~O)Rc?@SW`8A_W()OHMA^UAPHW`&l_Jj zI%2C)zKmVCtYwpqQk+of^o<;RI8P8w4ea1DsJM3sMk3t|gR|ZUO8*k61{?Z$X#rSp z8cTopbiv8(1c?svI5xL5nhh7EzG^EbgbII>wn+;CDRCt1hLk&rxZ&J!pxp+Die)Oq@X zcz*n6BgV@D31{}<#}B0fy6xf=@H`rYD(v(k0u78skp3Q9t1+g+b!%cvE4@l;_6;+s zRj;|zO}{sD;cp#_oIPV}>(00q#o0aBRM1v=_-_|`=`B_CVX!Bf<8W!8E`tc@Q>e|G zwI&z=FTw-*14%aX+${b)LmZ+EGY6R*A?y-`O#!2hQA0Ctu6#}OqygiIU_uZ=6|0<@ zZi~^lE+#Ar?~uDGp6n>j&$LNHNI?Phw8!ia2GCw%_15B~%ueuv83`NgcDL|SHJJ%$ ztv()2!zT6nAycGw zxM^&g1R-D{;C&k|QJl7iX)kRAs8|TJ4MJwnqOpVHBMvt;u||T>p^MH9DT&=E8trKNWQ#EZVKJ$4pVD?o zj8^3*7ui~CdOpys57M}G;7bv8uWt0+?13nP?e>w1yJeGA6$=$|iB7Ddm3lF+B17m} zVeX8Zz&8$G@ID?xtH#_CnxS_~-QwTl?k3#%Gy1Ggihg~)pqVyS9ntb!4VUN{2a#5^ zh=I=!G|z@nU?WUA{@Sbvp%Jzg%R{bvMq21_Y~(D=kz!+{M{l+2)s}NGw+pf|rrsm4 z)Vs;0>O&v`UyOrA?^h5`MyPMhjZPI|qP% zCSJ&gZP&o0waS#LaIX|=Qvtla(*@e!N>yS|4h%Er4s|189ua90YjN?A$iklY&6lFG zR+<2f&^s(p+s5%1#^r58x=20o?dhmN$zJP1t0K{!%C$zlmkCooXMyDtD*NdlN?BV7 z?xQN{%uQE7HvOoX;ZL@HNL$p4S8IbbM76BV0Z0v-4*swc3z00bvB&VB9!ZCHKxTH& z2-IZ`Owxr=PbTooTmob5Ly5t{JjAGjx= zc)TArqh1PxEHQ5VB~zxMqb-322shCSX?gErN3gO}BDO%^PTVB%OT-3~<-4A~GRZS` zq(fSW&Ss3LHCmDm81sY^rgFv1DgQqnwZfUCoiH{D7fj`{l~eU<9ccYsQmsH!R9ir1 z=oN`Q%}LaN-C2(i*}BW(2J5I2{RQ56X5gg^Y5@c(b)0H*xaA3ignL3-tS)&DF}89y zNLnLr+dxH?;ZJPFX^d^6Lvz0f^IdgI&)Oe7TATULzK0-!s7mvoPQ02z`wl(B$}ej) z+AOZX2p{&QRQ6SwQ^)|q&J;*P);3kGC-~~aGE-HSEO@l~VN8`>5rzE*g(b34c0XkaJgW`2cXjYYIz%IH;G3n$H5dTpNDmV-MdOnA?H2sMOv;CCY^fC zEMUM?TAsIrw83oT`E>RZYd-dYe@KFp>U^-&d=QmXHy8%;ysU<(H~_MNI!;qggF(&2 zAf`(whYxZQBR^TIZLtd)sp5RKD(~~v87ly=LnAHezBIbKi#8IJc|fE4n$ta}-xh1O z8B*tNtsS79v{R!Dr8^MPNiFreau7U+)l|4Q$qLYm`QAbPZ4BfDAIg+_;N!r9j}%6T z-*3Od?Q^%$Bb*N0QcUG2&a$W4-gA&Y@W7b&_uUbHv&J#q+gsSPh(G#)>CTmKhoT`p3%^m?@zM9mV=vYi+=JMMw3&o zpMPicwaC-PO;jPdA)7hHg{MwMVdq&M@Ea;h;p>14LRW^>^ot296yzfuaY6iuAMo-o zClkI;h*dm)l0~UC+h$qqP67hk=uQw4*FMtm2#R1ZOG~hz2hd>a*coq6+btm!&{Ls* zQE2zY89scRW;vc}n)bsIKv5^4$RC~z${?!Vo_QCgA1oSY$|iHH*`6e^1xUu8XIbd& z#~^K;(dzHX(JA#8Z`%daI^&kY^rn~1#mm&L#q6LTe7{my#;Sy<}FcpU(mQvztQ1dB%HIdu@!zI&K;Ue-@D5vujQ;L{6-soyR#qF#Z7iOGAn`ZPWQ)+hTa=? zHoV4Fj2jk@Y_01l5Fcm*)8^8XrqOt29KRhwmseeVqR5vb!@#!#Elthdie((7>l(#Z z!7~wuPR0MCCIhj6!Laj(aFqg-8d)Ej!zGX0*EP#^Pf@@mjc8xLQ*DyL@sqRfO%gVS zeLfp}pG;Wh@hvwl2Jp+1I@X}vPNfDOTb}9xWtB$|Dykf`#~e_FaUATiu_r23azMBH z9rp|hBB$pf_lx86;c~13SG2kx;QD&Z!q;~eV9o|c+pLYkkQ5KM1X*O+v#l0s~y&*I=RR~cF(O5!=OG+Nktvwy2)ordt@ zuvAgI>v0hy!#&F|D;F0CN=Jcw#L*J9_j09+JPoC8AJPca-tgLGCkksoW}OkguS!Xa zF80PL}y<#R71hWl6Ft?$LFuD^|Z4}p`9dI+fvw39x@fzLR zxfioh)XlpAPX5K@1&@afYeM#_9Y1}~Gf;d~yZ-ew^q2O_K)X8HQAg)ROZG^1QsSPd z&5nhL@&-h#eGD_>V;jrYc`i=lftifkIT80{1M&uH&d+g+$qu4~j4OG>g}%G^)qTfQ z<}T6He(cJ}ZQ}G|XjB%3O&l%nyr=ySqgR+XpwrM>4Ld+4KS&`^Eh^pd%=%8!!iu+} zTWO;wt>lv>Y#8H1Zd0j&wVzBk zan@)0Et4P=?8A~l4ZZnzGA^j9zk1{!F8oFp&6evkyld4Rm>k-x{V$j!Gc+^vsdseK zNGiHVMA&-5<{39bU6E`h(7TxnE_fJvBGz2F&=V6?lOrcn1OrB_DM9NU4@(j$Xl>`j*pZ{|rOQuq zq$^A3MYjCkMmw#pw&%pzd2Q6l)&wd)2`T|n@SjdaO{*G>a_Fm)sX(@aVuMy-^=0Yd zc7~x0cF~GUA#!`1To>OrFF1z|x?-J{{njpCKW4jT4|T9U9NRn-8eevgt3NlbU^2P$;!qLAo(?TftZ;QE)862;Db@Lwi$}XB zRR<1rs8vlX&jRHnA4McM{2+N90BvaLkHgjo9cA64g&mlXq-`QI0a|yzoP}Q}&VLKjimh zZfeGH$U@I|o%iR00pq*NA~JK&srZp6^jVJ~|8ToMOcI;;T)rZn;vGF=?7PUwUtMY! z;-@ru(@>XmR)oDlp&*8lX7z+^zYA_~uYOSFvqDXt7&15cej?mTJKbDB?*JvXU7&HJ#E``=tI*M&&w~FY)zO`r&TA)e!n3g-HQiC3ZkR*Y zSH~+)0&Ug0Qw{F;$P3IfA zFlnt^$j?fMlIN&sH$TEej7@YiGurk+A32-ez(O*q$^OIg@6TEemtHsmA8*(so8y(pi zT~PQaXE^{>jMb>o7YF*>pnAJDN=Pr<->^HbAYg<3dhonZk4m-0^-XFKs;2Lzr_&i* zCY_R{gkU(@U@pvppn%-@f=mVa$iBS!ToPw}>2tudBB$qCC+F~GJPs|P*OlO2zy7vj zbA%AP4KB|pmLdLc_TMEh>e2bnUdHDWjBMT;f{${<^!`^|3|YT&v(_brZ%eh6`iGX~ zCD?X$S*%pGMUCG-Kx178L<($5*uaQRfoI9gs0^?*o z8hE$>e)0Rg0|4o&o`PIOtcL5J9VNc8J>d(6qGf_K;dD7Oh}ByJ3s*S!mtsp0Wh5xa z@b29mXtP*w{#jARb}?GvsNWX&W7cqLSW|pVzU!O5E1djjxL5Ga0-SOnTun~%&YbD9 z=PUkYq9Cs{pSg2pE4O02pMIt`LeUT9#rzwx5N;LV^rWlrwJX8l*oB?GN@ zrUH~!Q|J4Cd5%Q%XCTK0BZ$SZu86A|Y&HD`7n|S=HuGGNVd%a^ddZDI!)M`PC_SS3rmokM9UlybYGeM+yh2 zun6s!s(5cGQe_rxjsxQM+Vafm!TS)gdwL1q;=L@|o8!-)c7YF5Ej@QC!glx{@qccg zU_T~@-m{Jk_j<`pWRle4EK=7f=7DdR$F^hRV@xYTTMa{YZc+EmeSR!_YvwY>?Afpz z%dr`k2G<<$XAXTX@fd5M)M4SA0&e|)m0sh+dW=`EuMGee|8T07X{>1Fsk|y3*E0ma z4jDd#>VfKn&C?Gz1F^)uFG>3q9q>n*IT$LKLb2R9EfWvpP%u*G)d)Mf;w<>O@Tf*} zr2rQ3Z8&LAs9Q~|IddSstiTwfM(0-F5qI(X~32gE?=3TWSyu2oa$WT3&>K*IQ^K1ES#B zAh^>RlU;bfN6oPk0$sSQBA+O3p0#ROWz5hrGq7dd@)`T(`5AwG9x!`{yn8k@LPYNU z(*K`^183=0rmugs6nFnzA;|w)ZDBiGXH#2eeH9OTQzza3R8!<6%*YKfAO!1=Jfe=! z;UEY}dhVY=@Fs(6aBdAp_OxExSnT0-#l!;;62>-t-b|0%u8H3b11h%$aCy|IDG_8a zXkc-_xBA%L7LUWQ_(KQh@EE2V?X+zv57SICp2*2fKj>XqlRd|)hD(?>K~DgqXNchr zOfwYrYS!4cKn_c6+rPy7e1>zcbECnhd9>UXzof*mNIf{y97zvUzR;s16g7dh`+o=w zwPDRD36eUSk;${2RsnfKQfz&l1s`(vgt3^V)a-FfH$4V@9%#BWc+<{hU;F=B+6WDj zl{U&+J6!Yx4AbvNMg1LT00^qxKTT4WZsq$`{{{S?VyPj3{J#AY$>^Vb2KryC?`Uf6 zL{ImhYf0kNe?$`X8_D!xwm^p{=-eS7afJwhPv;6U9--f$o9`O~RD__VF|fd$=@D1g zv(wK;mCtU3)EjKH5Uf`(VAqp&Ket6oBXA_#OM=ZM#o-8Aw_Sa;u2|0ENlzDreb5|K zkm?qiwlom$fK2cqjU(iwddXlBQzEt)Ro&7v_{ao&l5yCQJ`e5L@?b&}5tnpI*l$;m z52^RAC-c;OCw%+`xo5$;k6xYyW7Wf6jMy$yXS217!GC~#y=$*H;V)M$IN&U1WT)cb z(W~1u{gky@b(50_8D67=6G_phf0EdROtKLEJn650R& literal 0 HcmV?d00001 diff --git a/Common_Functions/atmo_p.m b/Common_Functions/atmo_p.m new file mode 100644 index 0000000..e7d33b0 --- /dev/null +++ b/Common_Functions/atmo_p.m @@ -0,0 +1,61 @@ +function P = atmo_p(alt, T, sum_n) +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +% Program: Atmospheric Pressure Calculation +% Author: Brent Lewis(RocketLion@gmail.com) +% University of Colorado-Boulder +% History: Original-1/10/2007 +% Revision-1/12/2007-Corrected for changes in Matlab versions +% for backward compatability-Many thanks to Rich +% Rieber(rrieber@gmail.com) +% Input: alt: Geometric altitude vector of desired pressure data[km] +% T: Temperature vector at given altitude points +% Required only for altitudes greater than 86 km[K] +% sum_n: Total number density of atmospheric gases[1/m^3] +% Output: P: Pressure vector[Pa] +% Note: Must compute altitudes below 86 km and above 86 km on two +% different runs to allow line up of altitudes and +% temperatures +% Examples: atmo_p(0) = 101325 Pa +% atmo_p(0:10) = Pressures between 0-10 km at 1 km increments +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +if nargin == 1 + T = []; + sum_n = []; +end + +% Constants +N_A = 6.022169e26; +g_0 = 9.80665; +M_0 = 28.9644; +R = 8.31432e3; +r_E = 6.356766e3; +% Geopotential/Geometric Altitudes used for Geometric Altitudes < 86 km +H = [0 11 20 32 47 51 71 84.852]; +Z = r_E*H./(r_E-H); +Z(8) = 86; + +% Defined temperatures/lapse rates/pressures/density at each layer +T_M_B = [288.15 216.65 216.65 228.65 270.65 270.65 214.65]; +L = [-6.5 0 1 2.8 0 -2.8 -2]/1e3; +P_ref = [1.01325e5 2.2632e4 5.4748e3 8.6801e2 1.1090e2 6.6938e1 3.9564]; + +% Preallocation of Memory +P = zeros(size(alt)); + +for i = 1 : length(alt) + Z_i = alt(i); + if isempty(sum_n) + index = find(Z>=Z_i)-1+double(Z_i==0); + index = index(1); + Z_H = r_E*Z_i/(r_E+Z_i); + if L(index) == 0 + P(i) = P_ref(index)*exp(-g_0*M_0*(Z_H-H(index))*1e3/(R*T_M_B(index))); + else + P(i) = P_ref(index)*(T_M_B(index)/... + (T_M_B(index)+L(index)*(Z_H-H(index))*1e3))^... + (g_0*M_0/(R*L(index))); + end + else + P(i) = sum_n(i)*R*T(i)/N_A; + end +end \ No newline at end of file diff --git a/Common_Functions/color_line.m b/Common_Functions/color_line.m new file mode 100644 index 0000000..249fc55 --- /dev/null +++ b/Common_Functions/color_line.m @@ -0,0 +1,43 @@ +function h = color_line(x, y, c, varargin) +% color_line plots a 2-D "line" with c-data as color +% +% h = color_line(x, y, c) +% by default: 'LineStyle','-' and 'Marker','none' +% +% or +% h = color_line(x, y, c, mark) +% or +% h = color_line(x, y, c, 'Property','value'...) +% with valid 'Property','value' pairs for a surface object +% +% in: x x-data +% y y-data +% c 3rd dimension for colouring +% mark for scatter plots with no connecting line +% +% out: h handle of the surface object + +% (c) Pekka Kumpulainen +% www.tut.fi + + +h = surface(... + 'XData',[x(:) x(:)],... + 'YData',[y(:) y(:)],... + 'ZData',zeros(length(x(:)),2),... + 'CData',[c(:) c(:)],... + 'FaceColor','none',... + 'EdgeColor','flat',... + 'Marker','none'); + +if nargin ==4 + switch varargin{1} + case {'+' 'o' '*' '.' 'x' 'square' 'diamond' 'v' '^' '>' '<' 'pentagram' 'p' 'hexagram' 'h'} + set(h,'LineStyle','none','Marker',varargin{1}) + otherwise + error(['Invalid marker: ' varargin{1}]) + end + +elseif nargin > 4 + set(h,varargin{:}) +end diff --git a/Common_Functions/coordinates_to_dxdy.m b/Common_Functions/coordinates_to_dxdy.m new file mode 100644 index 0000000..bf6ea6d --- /dev/null +++ b/Common_Functions/coordinates_to_dxdy.m @@ -0,0 +1,26 @@ +function [ sx, sy ] = coordinates_to_dxdy( long, lat ) +%constants +r_earth = 6371000; %(m) earth radius + +%calc +sx = 0; +sy = 0; +for n = 1:1:length(long) + if n==1 + else + dlong = long(n) - long(n-1); %(degrees) change in long position during this time slice + dlat = lat(n) - lat(n-1); %(degrees) change in lat position during this time slice + + meters_per_deg_long = (2*pi/360) * r_earth * cosd(lat(n)); %(m/deg) long + meters_per_deg_lat = (2*pi/360) * r_earth; %(m/deg) lat + + dx = meters_per_deg_long * dlong; %(m) + dy = meters_per_deg_lat * dlat; %(m) + + sx(n) = sx(end) + dx; %m + sy(n) = sy(end) + dy; %m + end +end + +end + diff --git a/Common_Functions/dxdy_to_coordinates.m b/Common_Functions/dxdy_to_coordinates.m new file mode 100644 index 0000000..99efce3 --- /dev/null +++ b/Common_Functions/dxdy_to_coordinates.m @@ -0,0 +1,26 @@ +function [ long, lat ] = dxdy_to_coordinates( sx, sy, long0, lat0 ) +%constants +r_earth = 6371000; %(m) earth radius + +%calc +long = long0; +lat = lat0; +for n = 1:1:length(sx) + if n==1 + else + dx = sx(n) - sx(n-1); %(m) change in x position during this time slice + dy = sy(n) - sy(n-1); %(m) change in y position during this time slice + + meters_per_deg_long = (2*pi/360) * r_earth * cosd(lat(end)); %(m/deg) long + meters_per_deg_lat = (2*pi/360) * r_earth; %(m/deg) lat + + dlong = dx / meters_per_deg_long; %(deg) + dlat = dy / meters_per_deg_lat; %(deg) + + long(n) = long(end) + dlong; + lat(n) = lat(end) + dlat; + end +end + +end + diff --git a/Common_Functions/plotxx.m b/Common_Functions/plotxx.m new file mode 100644 index 0000000..968f5e5 --- /dev/null +++ b/Common_Functions/plotxx.m @@ -0,0 +1,87 @@ +function [ax,hl1,hl2] = plotxx(x1,y1,x2,y2,xlabels,ylabels); +%PLOTXX - Create graphs with x axes on both top and bottom +% +%Similar to PLOTYY, but ... +%the independent variable is on the y-axis, +%and both dependent variables are on the x-axis. +% +%Syntax: [ax,hl1,hl2] = plotxx(x1,y1,x2,y2,xlabels,ylabels); +% +%Inputs: X1,Y1 are the data for the first line (black) +% X2,Y2 are the data for the second line (red) +% XLABELS is a cell array containing the two x-labels +% YLABELS is a cell array containing the two y-labels +% +%The optional output handle graphics objects AX,HL1,HL2 +%allow the user to easily change the properties of the plot. +% +%Example: Plot temperature T and salinity S +% as a function of depth D in the ocean +% +%D = linspace(-100,0,50); +%S = linspace(34,32,50); +%T = 10*exp(D/40); +%xlabels{1} = 'Temperature (C)'; +%xlabels{2} = 'Salinity'; +%ylabels{1} = 'Depth(m)'; +%ylabels{2} = 'Depth(m)'; +%[ax,hlT,hlS] = plotxx(T,D,S,D,xlabels,ylabels); +% +% +%The code is inspired from page 10-26 (Multiaxis axes) +%of the manual USING MATLAB GRAPHICS, version 5. +% +%Tested with Matlab 5.3.1 and above on PCWIN +% +%Author: Denis Gilbert, Ph.D., physical oceanography +%Maurice Lamontagne Institute, Dept. of Fisheries and Oceans Canada +%email: gilbertd@dfo-mpo.gc.ca Web: http://www.qc.dfo-mpo.gc.ca/iml/ +%November 1997; Last revision: 01-Nov-2001 +% +% +% Taken from MATLAB Central: +% +% http://www.mathworks.com/matlabcentral/fileexchange/317-plotxx-m + +if nargin < 4 + error('Not enough input arguments') +elseif nargin==4 + %Use empty strings for the xlabels + xlabels{1}=' '; xlabels{2}=' '; ylabels{1}=' '; ylabels{2}=' '; +elseif nargin==5 + %Use empty strings for the ylabel + ylabels{1}=' '; ylabels{2}=' '; +elseif nargin > 6 + error('Too many input arguments') +end + +if length(ylabels) == 1 + ylabels{2} = ' '; +end + +if ~iscellstr(xlabels) + error('Input xlabels must be a cell array') +elseif ~iscellstr(ylabels) + error('Input ylabels must be a cell array') +end + +hl1=line(x1,y1,'Color','k'); +ax(1)=gca; +set(ax(1),'Position',[0.12 0.12 0.75 0.70]) +set(ax(1),'XColor','k','YColor','k'); + +ax(2)=axes('Position',get(ax(1),'Position'),... + 'XAxisLocation','top',... + 'YAxisLocation','right',... + 'Color','none',... + 'XColor','r','YColor','k'); + +set(ax,'box','off') + +hl2=line(x2,y2,'Color','r','Parent',ax(2)); + +%label the two x-axes +set(get(ax(1),'xlabel'),'string',xlabels{1}) +set(get(ax(2),'xlabel'),'string',xlabels{2}) +set(get(ax(1),'ylabel'),'string',ylabels{1}) +set(get(ax(2),'ylabel'),'string',ylabels{2}) diff --git a/Common_Functions/savex.m b/Common_Functions/savex.m new file mode 100644 index 0000000..c8a199e --- /dev/null +++ b/Common_Functions/savex.m @@ -0,0 +1,143 @@ +function savex(varargin) +% % SAVEX +% % +% % save all variables that are present in the workspace EXCEPT what you +% % specify explicitly (by including the variable names or by using regular +% % expressions). +% % +% % Author : J.H. Spaaks +% % Date : April 2009 +% % MATLAB version : R2006a on Windows NT 6.0 32-bit +% % +% % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % % +% % +% % test1a = 0; +% % test2a = 2; +% % test3a = 4; +% % test4a = 6; +% % test5a = 8; +% % +% % test1b = 1; +% % test2b = 3; +% % test3b = 5; +% % test4b = 7; +% % test5b = 9; +% % +% % % This example saves all variables that are present in the workspace except +% % % 'test2a' and 'test5b'. 'test3' is ignored since there is no variable by +% % % that name: +% % savex('save-by-varname.mat','test2a','test3','test5b') +% % +% % % This example saves all variables that are present in the workspace except +% % % 'test4a', 'test4b' and 'test2b': +% % savex('save-by-regexp.mat','-regexp','test4[ab]','t[aeiou]st2[b-z]') +% % +% % % This example saves all variables that are present in the workspace except +% % % those formatted as Octave system variables, such as '__nargin__': +% % savex('no-octave-system-vars.mat','-regexp','^__+.+__$') +% % +% % % This example saves all variables that are present in the workspace except +% % % those which are specified using regular expressions, saving in ascii +% % % format. Supported options are the same as for SAVE. +% % savex('save-with-options.txt','-regexp','test4[ab]',... +% % 't[aeiou]st2[b-z]','-ascii') +% % +% % +% % +% % % clear +% % % load('save-by-varname.mat') +% % % +% % % clear +% % % load('save-by-regexp.mat') +% % % +% % % clear +% % % load('no-octave-system-vars.mat') +% % % +% % % clear +% % % load('save-with-options.txt','-ascii') +% % % + + +varList = evalin('caller','who'); +saveVarList = {}; + +if ismember(nargin,[0,1]) + % save all variables + saveVarList = varList + for u = 1:numel(saveVarList) + eval([saveVarList{u},' = evalin(',char(39),'caller',char(39),',',char(39),saveVarList{u},char(39),');']) + end + save('matlab.mat',varList{:}) + +elseif strcmp(varargin{2},'-regexp') + % save everything except the variables that match the regular expression + + optionsList = {}; + inputVars ={}; + for k=3:numel(varargin) + if strcmp(varargin{k}(1),'-') + optionsList{1,end+1} = varargin{k}; + else + inputVars{1,end+1} = varargin{k}; + end + end + + + for k=1:numel(varList) + + matchCell = regexp(varList{k},inputVars,'ONCE'); + + matchBool = repmat(false,size(matchCell)); + for m=1:numel(matchCell) + if ~isempty(matchCell{m}) + matchBool(m) = true; + end + end + if ~any(matchBool) + saveVarList{end+1} = varList{k}; + end + end + + for u = 1:numel(saveVarList) + eval([saveVarList{u},' = evalin(',char(39),'caller',char(39),',',char(39),saveVarList{u},char(39),');']) + end + + save(varargin{1},saveVarList{:},optionsList{:}) + + + +elseif ischar(varargin{1}) + % save everything except the variables that the user defined in + % varargin{2:end} + optionsList = {}; + inputVars = {}; + for k=2:numel(varargin) + if strcmp(varargin{k}(1),'-') + optionsList{1,end+1} = varargin{k}; + else + inputVars{1,end+1} = varargin{k}; + end + end + + for k=1:numel(varList) + + if ~ismember(varList{k},inputVars) + + saveVarList{end+1} = varList{k}; + + end + + end + + for u = 1:numel(saveVarList) + eval([saveVarList{u},' = evalin(',char(39),'caller',char(39),',',char(39),saveVarList{u},char(39),');']) + end + + save(varargin{1},saveVarList{:},optionsList{:}) + +else + error('Unknown function usage.') +end + + + diff --git a/Common_Functions/sphere_CD.m b/Common_Functions/sphere_CD.m new file mode 100644 index 0000000..9fd4cc5 --- /dev/null +++ b/Common_Functions/sphere_CD.m @@ -0,0 +1,8 @@ +function CD = sphere_CD( RE ) +% references: +% http://www.chem.mtu.edu/~fmorriso/DataCorrelationForSphereDrag2013.pdf + +CD = (24./RE) + (2.6*(RE./5))/(1+(RE./5).^1.52) + 0.411*(RE/263000).^-7.94./(1+(RE/263000).^-8) + RE.^0.8/461000; + +end + diff --git a/Common_Functions/stdatmo.m b/Common_Functions/stdatmo.m new file mode 100644 index 0000000..25e4ea1 --- /dev/null +++ b/Common_Functions/stdatmo.m @@ -0,0 +1,405 @@ +function [rho,a,temp,press,kvisc,ZorH]=stdatmo(H_in,Toffset,Units,GeomFlag) +% STDATMO Find gas properties in earth's atmosphere. +% [rho,a,T,P,nu,ZorH]=STDATMO(H,dT,Units,GeomFlag) +% +% STDATMO by itself gives the atmospheric properties at sea level on a +% standard day. +% +% STDATMO(H) returns the properties of the 1976 Standard Atmosphere at +% geopotential altitude H (meters), where H is a scalar, vector, matrix, +% or ND array. +% +% STDATMO(H,dT) returns properties when the temperature is dT degrees +% offset from standand conditions. H and dT must be the same size or else +% one must be a scalar. +% +% STDATMO(H,dT,Units) specifies units for the inputs outputs. Options are +% SI (default) or US (a.k.a. Imperial, English). For SI, set Units to [] +% or 'SI'. For US, set Units to 'US'. Input and output units may be +% different by passing a cell array of the form {Units_in Units_out}, +% e.g. {'US' 'SI'}. Keep in mind that dT is an offset, so when converting +% between Celsius and Fahrenheit, use only the scaling factor (dC/dF = +% dK/dR = 5/9). Units are as follows: +% Input: SI (default) US +% H: Altitude m ft +% dT: Temp. offset °C/°K °F/°R +% Output: +% rho: Density kg/m^3 slug/ft^3 +% a: Speed of sound m/s ft/s +% T: Temperature °K °R +% P: Pressure Pa lbf/ft^2 +% nu: Kinem. viscosity m^2/s ft^2/s +% ZorH: Height or altitude m ft +% +% STDATMO(H,dT,u), where u is a structure created by the UNITS function, +% accepts variables of the DimensionedVariable class as inputs. Outputs +% are of the DimensionedVariable class. If a DimensionedVariable is not +% provided for an input, STDATMO assumes SI input. +% +% STDATMO(H,dT,Units,GeomFlag) with logical input GeomFlag returns +% properties at geometric altitude input H instead of the normal +% geopotential altitude. +% +% [rho,a,T,P,nu]=STDATMO(H,dT,...) returns atmospheric properties the +% same size as H and/or dT (P does not vary with temperature offset and +% is always the size of H) +% +% [rho,a,T,P,nu,ZorH]=STDATMO(H,...) returns either geometric height, Z, +% (GeomFlag not set) or geopotential height, H, (GeomFlag set). +% +% Example 1: Find atmospheric properties at every 100 m of geometric +% height for an off-standard atmosphere with temperature offset varying +% +/- 25°C sinusoidally with a period of 4 km. +% Z = 0:100:86000; +% [rho,a,T,P,nu,H]=stdatmo(Z,25*sin(pi*Z/2000),'',true); +% semilogx(rho/stdatmo,H/1000) +% title('Density variation with sinusoidal off-standard atmosphere') +% xlabel('\sigma'); ylabel('Altitude (km)') +% +% Example 2: Create tables of atmospheric properties up to 30000 ft for a +% cold (-15°C), standard, and hot (+15°C) day with columns +% [h(ft) Z(ft) rho(slug/ft3) sigma a(ft/s) T(R) P(psf) µ(slug/ft-s) nu(ft2/s)] +% using 3-dimensional array inputs. +% [~,h,dT]=meshgrid(0,-5000:1000:30000,-15:15:15); +% [rho,a,T,P,nu,Z]=stdatmo(h,dT*9/5,'US',0); +% Table = [h Z rho rho/stdatmo(0,0,'US') T P nu.*rho nu]; +% format short e +% ColdTable = Table(:,:,1) +% StandardTable = Table(:,:,2) +% HotTable = Table(:,:,3) +% +% Example 3: Use the unit consistency enforced by the DimensionedVariable +% class to find the SI dynamic pressure, Mach number, Reynolds number, and +% stagnation temperature of an aircraft flying at flight level FL500 +% (50000 ft) with speed 500 knots and characteristic length of 80 inches. +% u = units; +% V = 500*u.kts; c = 80*u.in; +% [rho,a,T,P,nu]=stdatmo(50*u.kft,[],u); +% Dyn_Press = 1/2*rho*V^2; +% M = V/a; +% Re = V*c/nu; +% T0 = T*(1+(1.4-1)/2*M^2); +% +% This atmospheric model is not recommended for use at altitudes above +% 86 km geometric height (84852 m/278386 ft geopotential) and returns NaN +% for altitudes above 90 km geopotential. +% +% See also DENSITYALT, ATMOSISA, ATMOSNONSTD, DENSITYALT, +% UNITS, DIMENSIONEDVARIABLE. +% http://www.mathworks.com/matlabcentral/fileexchange/39325 +% http://www.mathworks.com/matlabcentral/fileexchange/38977 +% +% [rho,a,T,P,nu,ZorH]=STDATMO(H,dT,Units,GeomFlag) + +% About: +%{ +Author: Sky Sartorius + www.mathworks.com/matlabcentral/fileexchange/authors/101715 + +References: ESDU 77022 + www.pdas.com/atmos.html +%} + +if nargin >= 3 && isstruct(Units) + U = true; + u = Units; + if isa(H_in,'DimensionedVariable') + H_in = H_in/u.m; + end + if isa(Toffset,'DimensionedVariable') + Toffset = Toffset/u.K; + end + + Units = 'si'; +else + U = false; +end + +if nargin == 0 + H_in = 0; +end + +if nargin < 2 || isempty(Toffset) + Toffset = 0; +end + +if nargin <= 2 && all(H_in(:) <= 11000) %quick troposphere-only code + TonTi=1-2.255769564462953e-005*H_in; + press=101325*TonTi.^(5.255879812716677); + temp = TonTi*288.15 + Toffset; + rho = press./temp/287.05287; + + if nargout > 1 + a = sqrt(401.874018 * temp); + if nargout >= 5 + kvisc = (1.458e-6 * temp.^1.5 ./ (temp + 110.4)) ./ rho; + if nargout == 6 % Assume Geop in, find Z + ZorH = 6356766*H_in./(6356766-H_in); + end + end + end + return +end + +% index Lapse rate Base Temp Base Geopo Alt Base Pressure +% i Ki (°C/m) Ti (°K) Hi (m) P (Pa) +D =[1 -.0065 288.15 0 101325 + 2 0 216.65 11000 22632.0400950078 + 3 .001 216.65 20000 5474.87742428105 + 4 .0028 228.65 32000 868.015776620216 + 5 0 270.65 47000 110.90577336731 + 6 -.0028 270.65 51000 66.9385281211797 + 7 -.002 214.65 71000 3.9563921603966 + 8 0 186.94590831019 84852.0458449057 0.373377173762337 ]; + +% Constants +R=287.05287; %N-m/kg-K; value from ESDU 77022 +% R=287.0531; %N-m/kg-K; value used by MATLAB aerospace toolbox ATMOSISA +gamma=1.4; +g0=9.80665; %m/sec^2 +RE=6356766; %Radius of the Earth, m +Bs = 1.458e-6; %N-s/m2 K1/2 +S = 110.4; %K + +K=D(:,2); %°K/m +T=D(:,3); %°K +H=D(:,4); %m +P=D(:,5); %Pa + +temp=zeros(size(H_in)); +press=temp; +hmax = 90000; + +if nargin < 3 || isempty(Units) + Uin = false; + Uout = Uin; +elseif isnumeric(Units) || islogical(Units) + Uin = Units; + Uout = Uin; +else + if ischar(Units) %input and output units the same + Unitsin = Units; Unitsout = Unitsin; + elseif iscell(Units) && length(Units) == 2 + Unitsin = Units{1}; Unitsout = Units{2}; + elseif iscell(Units) && length(Units) == 1 + Unitsin = Units{1}; Unitsout = Unitsin; + else + error('Incorrect Units definition. Units must be ''SI'', ''US'', or 2-element cell array') + end + + if strcmpi(Unitsin,'si') + Uin = false; + elseif strcmpi(Unitsin,'us') + Uin = true; + else error('Units must be ''SI'' or ''US''') + end + + if strcmpi(Unitsout,'si') + Uout = false; + elseif strcmpi(Unitsout,'us') + Uout = true; + else error('Units must be ''SI'' or ''US''') + end +end + +% Convert from imperial units, if necessary. +if Uin + H_in = H_in * 0.3048; + Toffset = Toffset * 5/9; +end + + + +% Convert from geometric altitude to geopotental altitude, if necessary. +if nargin < 4 + GeomFlag = false; +end +if GeomFlag + Hgeop=(RE*H_in)./(RE+H_in); +else + Hgeop=H_in; +end + +n1=(Hgeop<=H(2)); +n2=(Hgeop<=H(3) & Hgeop>H(2)); +n3=(Hgeop<=H(4) & Hgeop>H(3)); +n4=(Hgeop<=H(5) & Hgeop>H(4)); +n5=(Hgeop<=H(6) & Hgeop>H(5)); +n6=(Hgeop<=H(7) & Hgeop>H(6)); +n7=(Hgeop<=H(8) & Hgeop>H(7)); +n8=(Hgeop<=hmax & Hgeop>H(8)); +n9=(Hgeop>hmax); + +% Troposphere +if any(n1(:)) + i=1; + TonTi=1+K(i)*(Hgeop(n1)-H(i))/T(i); + temp(n1)=TonTi*T(i); + PonPi=TonTi.^(-g0/(K(i)*R)); + press(n1)=P(i)*PonPi; +end + +% Tropopause +if any(n2(:)) + i=2; + temp(n2)=T(i); + PonPi=exp(-g0*(Hgeop(n2)-H(i))/(T(i)*R)); + press(n2)=P(i)*PonPi; +end + +% Stratosphere 1 +if any(n3(:)) + i=3; + TonTi=1+K(i)*(Hgeop(n3)-H(i))/T(i); + temp(n3)=TonTi*T(i); + PonPi=TonTi.^(-g0/(K(i)*R)); + press(n3)=P(i)*PonPi; +end + +% Stratosphere 2 +if any(n4(:)) + i=4; + TonTi=1+K(i)*(Hgeop(n4)-H(i))/T(i); + temp(n4)=TonTi*T(i); + PonPi=TonTi.^(-g0/(K(i)*R)); + press(n4)=P(i)*PonPi; +end + +% Stratopause +if any(n5(:)) + i=5; + temp(n5)=T(i); + PonPi=exp(-g0*(Hgeop(n5)-H(i))/(T(i)*R)); + press(n5)=P(i)*PonPi; +end + +% Mesosphere 1 +if any(n6(:)) + i=6; + TonTi=1+K(i)*(Hgeop(n6)-H(i))/T(i); + temp(n6)=TonTi*T(i); + PonPi=TonTi.^(-g0/(K(i)*R)); + press(n6)=P(i)*PonPi; +end + +% Mesosphere 2 +if any(n7(:)) + i=7; + TonTi=1+K(i)*(Hgeop(n7)-H(i))/T(i); + temp(n7)=TonTi*T(i); + PonPi=TonTi.^(-g0/(K(i)*R)); + press(n7)=P(i)*PonPi; +end + +% Mesopause +if any(n8(:)) + i=8; + temp(n8)=T(i); + PonPi=exp(-g0*(Hgeop(n8)-H(i))/(T(i)*R)); + press(n8)=P(i)*PonPi; +end + +if any(n9(:)) + warning('One or more altitudes above upper limit.') + temp(n9)=NaN; + press(n9)=NaN; +end + +temp = temp + Toffset; + +rho = press./temp/R; + +if nargout >= 2 + a = sqrt(gamma * R * temp); + if nargout >= 5 + kvisc = (Bs * temp.^1.5 ./ (temp + S)) ./ rho; %m2/s + if nargout == 6 + if GeomFlag % Geometric in, ZorH is geopotential altitude (H) + ZorH = Hgeop; + else % Geop in, find Z + ZorH = RE*Hgeop./(RE-Hgeop); + end + end + end +end + +if Uout %convert to imperial units if output in imperial units + rho = rho / 515.3788; + if nargout >= 2 + a = a / 0.3048; + temp = temp * 1.8; + press = press / 47.88026; + if nargout >= 5 + kvisc = kvisc / 0.09290304; + if nargout == 6 + ZorH = ZorH / 0.3048; + end + end + end +end + +if U + rho = rho*u.kg/(u.m^3); + if nargout >= 2 + a = a*u.m/u.s; + temp = temp*u.K; + press = press*u.Pa; + if nargout >= 5 + kvisc = kvisc*u.m^2/u.s; + if nargout == 6 + ZorH = ZorH*u.m; + end + end + end +end + +end + +% Credit for elements of coding scheme: +% cobweb.ecn.purdue.edu/~andrisan/Courses/AAE490A_S2001/Exp1/ + +% Revision history: +%{ +V1.0 5 July 2010 +V1.1 8 July 2010 + Update to references and improved input handling +V2.0 12 July 2010 + Changed input ImperialFlag to Units. Units must now be a string or cell + array {Units_in Units_out}. Version 1 syntax works as before. + Two examples added to help +V2.1 15 July 2010 + Changed help formatting + Sped up code - no longer caclulates a or nu if outputs not specified. + Also used profiler to speed test against ATMOSISA, which is + consistently about 5 times slower than STDATMO + 17 July 2010 + Cleaned up Example 1 setup using meshgrid + 26 July 2010 + Switched to logical indexing, which sped up running Example 1 + significantly(running [rho,a,T,P,nu,h]=stdatmo(Z,dT,'US',1) 1000 times: + ~.67s before, ~.51s after) +V3.0 7 August 2010 + Consolodated some lines for succintness Changed Hgeop output to be + either geopotential altitude or geometric altitude, depending on which + was input. Updated help and examples accordingly. +V3.1 27 August 2010 + Added a very quick, troposhere-only section +V3.2 23 December 2010 + Minor changes, tested on R2010a, and sinusoidal example added +V4.0 6 July 2011 + Imperial temp offset now °F/°R instead of °C/°K +V4.1 12 Sep 2012 + Added ZorH output support for quick troposphere calculation + uploaded +V4.2 + tiny changes to help and input handling + nov 2012: some :s added to make use of any() better + added see alsos + uploaded +V5.0 + STDATMODIM wrapper created that takes DimensionedVariable input + uploaded 5 Dec 2012 +V6.0 + STDATMODIM functionality integrated into STDATMO; example three changed + for illustration. +%} \ No newline at end of file diff --git a/Common_Functions/stdatmo_H2.m b/Common_Functions/stdatmo_H2.m new file mode 100644 index 0000000..dbaccd4 --- /dev/null +++ b/Common_Functions/stdatmo_H2.m @@ -0,0 +1,22 @@ +function [rho_H2,SOS_H2,T_H2,P_H2]=stdatmo_H2(altitude) +%provides hydrogen properties at given altitude +%assumes the temperature and pressure of the air and H2 are the same +%all units are SI + +%% +%constants +R_air = 287.058; %specific gas constant air (SI) +R_H2 = 4124; %specific gas constant hydrogen (SI) +CP_H2 = 14.32; %kj/kg-K +CV_H2 = 10.16; %kj/kg-K +gamma_H2 = CP_H2 / CV_H2; %ratio of specific heats + +%standard atmosphere of the ambient air +[~,~,T_air,P_air]=stdatmo(altitude); %SI (standard atmosphere) +T_H2 = T_air; +P_H2 = P_air; +rho_H2 = P_H2 / (R_H2 * T_H2); +SOS_H2 = sqrt(gamma_H2 * R_H2 * T_H2); + +end + diff --git a/Common_Functions/transpose_arrayOfStructs.m b/Common_Functions/transpose_arrayOfStructs.m new file mode 100644 index 0000000..437cd86 --- /dev/null +++ b/Common_Functions/transpose_arrayOfStructs.m @@ -0,0 +1,16 @@ +function structOfArrays = transpose_arrayOfStructs( arrayOfStructs ) + +fields = fieldnames(arrayOfStructs(1)); + +for i = 1:length(fields) + + field = fields{i}; + + % When field = 'a', this is like saying: + % structOfArrays.a = [arrayOfStructs.a]; + structOfArrays.(field) = [arrayOfStructs.(field)]; + +end + +end + diff --git a/Common_Functions/volRevolve.m b/Common_Functions/volRevolve.m new file mode 100644 index 0000000..819983e --- /dev/null +++ b/Common_Functions/volRevolve.m @@ -0,0 +1,157 @@ +function vol = volRevolve(R,Z) +%VOLREVOLVE calculates the volume of a polygon revolved around the Z-axis. +% +% CALLING SEQUENCE / USAGE: +% +% vol = volRevolve(R,Z) +% +% INPUTS: +% +% R One-dimensional vector (n-by-1 or 1-by-n) of R points +% Z One-dimensional vector (n-by-1 or 1-by-n) of Z points +% +% OUTPUTS: +% +% vol Volume of solid formed by revolving polygon formed by [R,Z] +% lists of points around the Z-axis. Units are the square of the +% units of the inputs. +% +% DESCRIPTION and NOTES: +% +% - R and Z vectors must be in order, counter-clockwise around the area +% being defined. If not, this will give the volume of the +% counter-clockwise parts, minus the volume of the clockwise parts. +% +% - It does not matter if the curve is open or closed - if it is open +% (last point doesn't overlap first point), this function will +% automatically close it. +% +% - Based on Advanced Mathematics and Mechanics Applications with MATLAB, +% 3rd ed., by H.B. Wilson, L.H. Turcotte, and D. Halpern, +% Chapman & Hall / CRC Press, 2002, e-ISBN 978-1-4200-3544-5. +% See Chapter 5, Section 5.4, doi: 10.1201/9781420035445.ch5 +% +% Function originally created 2012-05-03 by Geoffrey M. Olynyk. +% (email: geoff at olynyk.name or golynyk at psfc.mit.edu) +% For details and notes, see notebook of G.M. Olynyk, dated 2012-05-03. + +% UPDATE HISTORY: +% +% 2012-05-03: Function originally created. -G.M. Olynyk + +% Check to make sure they're the same length, throw an error if not: +if ~(numel(R) == numel(Z)), + err = MException('volRevolve:UnequalLengths', ... + 'Input coordinate vectors (R,Z) do not have same length') ; + throw(err) ; +end + +% Now that we know they're the same length, figure out what that length is: +n = numel(R) ; + +% You have to have at least three points to define a polygon. +% Check that this is true and throw an error if it's not: +if n < 3, + err = MException('volRevolve:NotEnoughPoints', ... + 'Input coordinate vectors must have at least three points') ; + throw(err) ; +end + +% Cast inputs to doubles: +R = double(R) ; +Z = double(Z) ; + +% Check if last point overlaps first point; if it doesn't, add another +% point that does overlap the first point. (Note that we have to compare +% using a tolerance because these are floating-point values): +tol = 1.e-5 ; % (10 µm when using metres as input) + +if (abs(R(end) - R(1)) < tol) && (abs(Z(end) - Z(1)) < tol), + isClosed = true ; +else + isClosed = false ; +end + +% If it's not a closed curve, add a point to make it closed: +if ~isClosed, + R(n+1) = R(1) ; + Z(n+1) = Z(1) ; +end + +clear n ; +clear isClosed ; + +% Now, from the Wilson, Turcotte, and Halpern book, we note that if we have +% a closed curve defined by R(s), Z(s), 0 <= s <= 1, then the volume of the +% solid formed by revolving that curve around the Z axis is: +% +% V = (2pi/3) * int( R(s) * [R(s)Z'(s) - Z(s)R'(s)] ds, s = 0..1) +% +% where R'(s) = dR/ds; Z'(s) = dZ/ds. Since the boundary of our polygon is +% piecewise linear, these R'(s) and Z'(s) are constant on each line +% segment, and can be pre-calculated quickly. Note that the s values are +% equally spaced on each point. So, for example, if there are four line +% segments, the s values at points 1, 2, 3, 4, 5 are 0.00, 0.25, 0.50, +% 0.75, 1.00 respectively. (Note point 5 is on top of point 1.) + +Rp = R(1:end-1) ; +Zp = Z(1:end-1) ; + +dR = (R(2:end) - Rp) ; +dZ = (Z(2:end) - Zp) ; + +% We can do an analytic integral for each of nSeg line segments. It has +% four parts, we calculate each of these as v1, v2, v3, v4, and then add +% them together to get the total integral. (Everything is worked out in +% notebook of G.M. Olynyk dated 2012-05-03.) + +v1 = dR .* dZ .* Rp ; +v2 = 2 * dZ .* Rp.^2 ; +v3 = (-1) * dR.^2 .* Zp ; +v4 = (-2) * dR .* Rp .* Zp ; + +V = (pi/3) * (v1 + v2 + v3 + v4) ; +vol = sum(V) ; + + + +% Verification for circular toroid of major radius R0, minor radius a +% Volume is 2 * pi^2 * R0 * a^2. Run this code: + +% clear all +% R0 = 5 ; +% a = 1 ; +% npoints = 100 ; +% theta = 2*pi*[0:1:npoints-1]'/double(npoints-1) ; +% R = R0 + a*cos(theta) ; +% Z = a*sin(theta) ; +% vol_analytic = 2 * pi^2 * R0 * a^2 ; +% >> 98.6960 +% vol = volRevolve(R,Z) ; +% >> 98.6298 (6.7e-04 relative error) + +% Do it again with npoints = 1000, get: +% >> 98.6954 (6.6e-06 relative error) + +% As expected, it's always slightly small because the polygon inscribes the +% circle. + + + +% Verification for washer (rectangular toroid), with the radius of the +% 'hole' in the washer being a, and the outer radius of the washer being b. +% (Thus the width of the metal cross section is b-a.) The height of the +% washer is h. Then the volume is pi * (b^2 - a^2) * h. Run this code: +% +% clear all +% a = 1 ; +% b = 2 ; +% h = 10 ; +% R = [a; b; b; a; a] ; +% Z = [0; 0; h; h; 0] ; +% vol_analytic = pi * (b^2 - a^2) * h ; +% >> 94.2478 +% vol = volRevolve(R,Z) ; +% >> 94.2478 + +end % function \ No newline at end of file diff --git a/Common_Functions/wraptopi.m b/Common_Functions/wraptopi.m new file mode 100644 index 0000000..9c73243 --- /dev/null +++ b/Common_Functions/wraptopi.m @@ -0,0 +1,28 @@ +function a = wraptopi(a, a_center ) +% function a = wrap(a,a_center) +% +% Wraps angles to a range of 2*pi. +% Inverse of Matlab's "unwrap", and better than wrapToPi ( which has +% redundant [-pi,pi]) +% Optional input "a_center" defines the center angle. Default is 0, giving +% angles from (-pi,pi], chosen to match angle(complex(-1,0)). Maximum +% possible value is pi. + +% T.Hilmer, UH +% 2010.10.18 version 2 +% removed code from version 1. Have not bug-checked second input +% "a_center" + +if nargin < 2, a_center = 0; end + +% new way +a = mod(a,2*pi); % [0 2pi) + +% shift +j = a > pi - a_center; +a(j) = a(j) - 2*pi; +j = a < a_center - pi; +a(j) = a(j) + 2*pi; + + + diff --git a/GUI.fig b/GUI.fig new file mode 100644 index 0000000000000000000000000000000000000000..e2a42b9453dd78f75895932770fa32dc5020842c GIT binary patch literal 5729 zcma)AXEYlQ*SCsNtF+XNTD7TCJ5@^Uy?3qJVsC9xwbiV>OHtIGN$pMjWADZmD+!H6 z#M}S-e0aV-_uO;t_;AlX_nhBd+Db-RN-qJDB4Pk-B_m;HR}Uv4fS!lFzq7Zmrz}8E z*;HLe2%zliWbg0f2=Mlj1sDc+0d%|r0b()$aY@-{&t;!V1D=YBKLh;Fw822|A0kxO zrunOx{I!D%j?JgMSza*|OAl*Ps9VSA5!}UmNPkV|d$;Tkcatj-w*jh@mtVMPCnabE z&IvD^pR7~g`o~@LlY?jp{%(OZT92C=0NP3;eobCVZmb_O7|SKmEdgH-U|xkAl|I&8 z>5145@!q(e8VQqoL`Hx^-66Qsr#6^D%&AbmQ@Z45BTp4pZUMG>SaH4E*L;MYfc@@f zH{fOd)mor;{sck0%Oq|m)##wVlVGyPNw3UzATClQSaPCY>D=^- zgI+TmpQI0;aw#8$&7eO7bj8BTOltK>=in4WwAoUS@EO zuLNfQM4)8vrK4myi%f}S3pD*)8?1tR6t>C3^sB0hz`2wd;vtABj&CT+BwQC^%W$iW z*nDuSjEESzHNt&{ktJ{l;1{EkIKJY-GZP$u52{ka^kfHsho+D#ch+Yj%zVlD!}{Ab zpoI<^4VSApT+2C7f#QhRCtNX%ksFnp$++=t^mt5w!+zw)Cu?&Rk9_(D1IA@YOrQ$b z&(r94#av7eg1>aN4*7cl2`{UU9~KUJ6eKNvDh9X=)eY^8v7sK#D`D&ne<@011uGMZ zOpqo^@MdGRWgNW84Cx_xb{LTN{H2XSc>kyyddgfWYV#?0qx*V(V_jyTo^r)jd~(dh zv+1%({B^pbbGF(vydHK9sD8*-er1}HGOGHlgf{ALVW794v(73c<`8* z+7U*l zp3xHPSjQwXzKwFHZ*M@=suVau!(McqMKsD?Bg{uV(%5)|Q>Wgcn~C^#)&u#V(AhBY@{PBuy`PgpMIGlR5;RSE!&Md1^u3SCdn1uk z8H$WhJhf92%5|2wq+NFZgrVMG+G77%Rw1;mq3e6O-#1VtYf*WQj>_gW4zb@k9fCmN zDo$6EL?Tna; zz6ob@X~&6WGP#f&Jw1A);axU^)tEUmx>t2I?U9q!huKz=7(AO0!ABJI>yFpqEtbEy z)y*p^@fIry6K3!%?6yclzw+ijk*jS>tzEhdpRNg&(3uf{pEm5LuiBD>(!n(VV`gTn_4?Z)AYKeUPN#-Cp< zvc8@UY6y`$^~Qw>pn`Uf1&4!b$!JVz5{6TOg@3sprJ(&a5Fx=3P{yFoMD0Gb5pyXuZ_}`^(NBGm7B=-lyeC|33*dzK4&BPUK4!H zu{RY(NI-CVdu*<&ORvOG+@^H2*I>bEWE=%53w2cGg=|4yn_}-^9WMs$FKPeFi&86*iNI zb0$11?rTw%X-dn`fSop+k1<~PU1ob&wPt*9l-)#>+2+tE^Tkb*l2EInjR$M<#v+y`IS=1S_#z^ha8nZK62)>5JCLw4Plq^(Dv>6fut!;6XTC`TJv8awTV z*H&^p@EtWx>pv|xc^u)HcKzkf4Z6>%m5{31U`i}Z+6K=-*mSG*;t9ji5!#C|Q)mBsdEbD8lvZ`S+ zCLkd;&$0h)ng$}^RgZi8MPqLj}%WqZ7`$uFkPrx`bNY^ z7p?M&kVjSrj>!r+!ZNW4vkn+Ge;&!p76@8?aY=e-VBb?v7ibj8UVJ01XlYPaWiID^ zYso%1s{0{%l+{Xu${y=t)gsHZ`HZq^1Q zN`Vc5H3bZn&ob`Qmc!l$K?6|jBMH>`-;5uO zsS5aYNEd(iU}!65G!YPGpPb&cuhCiw=mKUqTu)+;xh{^m-i^QWE#t{aTJ{y^=p>n! zw?{mRD<%LGW_w^#*yqAS6E7pm@9;pQ&LWOi12L(H;f-E%ap+P_?`NL$44UfvqqLvL zE7Q$dSSQLr+K-oR`Z=7RO?iY1<=UTBwHQnJ@hl6Wn-3gjIwshMm8skddA(Q`!hjGz z-=p8YYSf+i^#MflE&6K$5_oFqAdcEL7taVdbwq}UoT@d*5toFy;@Z0!o;gi1O$!ru zMR=(G0(8@UyY#H~xqcXdQy>k!84PBJ*^m0kT?t_iPzh0%S&gdxNvxp-fyvk_SJovo zEI*AQ7nc+76Is-Db>@HkR^WHgXDCJHp9gI1t)X;Wwfzq`UrJ&(I=$H7M{=RepD11% zx9|!o@px1_E3O*xcGd9rhGnY9@nIenI+bo)m(IVeA;}%*dv>?`tWI!NSqat$r!KD; z(xkE&=6MYN_bV-Y!@2jBPMT@AzW*A?@l$Xl}#NU>l$UeA zssPbdcl%3za3F0z@Z9`jHVDxX@-;C(DsB>(k`E;CUX%BbXJxqdjrvpqu_fz+C}Anm z(W50d2o#@c>QIB>aU3-LEt0YQyV-Rro(n+DN@>6oNz@C1&2h|9 z8t4{;z9^7>+RSlz8hp5tukFa0`ivO}kFbT<&h$u#5E3|F3B9+=lzBXAmx%Mrn2}Kw zjj3g7<`hx)MTHB5<8c{G9JQOkQ0#hn6-$llryt+kYWYs)Ga$gM%jt@P-Y<=OKc{9G z7W!6x#9Ent^Z&8;OXSTcsJ1Wc<(b7tn`PaP>{pSIEE)M4p9P0@oIBTEej3s!nE<~d zfnWbhI$j!WbTC2c+SyO+e%N{6Qu7GA#}h#mW>)^t?@9b&{A!j74MKVJ=bUJC&>AY-owJ$QT$hhi1QZ<1lnSHO!R%bj9{(- z-#rOkH4P^gPLrv_%?}o4bUK-grV=lXPa(=e=da=k* zNM^Df0nkx>BQ%u%#=i<4+ut5vdkYeYyCT1b{6TovRry4^|0#_V^|)0=hzU0&GWQhD z?{dYz(O!ibdbP9~oAk3ynV%9I^m&Tv%^V?J8&jKJrm9_VHtz?YrT5v>K@~${U{O>z zFGkpByFB`Pr<$s;MX^sufl>h zA|GBfB)T;SGcU=xjf`|3#i#Di*~7q3hBUVTyBc4ReOC`=bC3m5hMI_^3qsea^38!H z16btN@3w{&|7w?L=w*ZNw}(5IS+mF=46se5qYz&X3CwGzm~GNtZ%=o7K}{KMDe|_H zox2jdOnJWP9-RHp$tuY5RgUR$JbYleuW`rhvt+hHO~Ab|xVFc+qO*Y9FJqe})}t`# zH;@x8>`Cs$DpKZh^(sVR3LIvDBZ+i{TsBgLOxJVpUOhi-21FHGCwe3e$>7=fr%A89)3eBZXIP^_=h5ccF1(C;Pe^28ri<-nQ4Ib+x|oj+nP9NPWBx z7kc%ds72+Tk016WwpwXz$q?&=F>*CbZOWcIVo6pTK4PoJYs3K=_OlsBkau_F&!*HQA{2PkQ>_?bH!fU%ZfyA?Znd!ZUJ3ZM z{P3FS+p^Ll4?^zk_O1OX7);5y+pZ&2ySuoPdDH#^+={=I564sf+qX_kCcL087+RrH zEhX8F(biJccN6qD$c;UYo*DDskNDJ5O{sY@Sj4Mb@MwY^>Nt?sPvdVm$k(_-jHO#Q zf)Ufgh|R>F7jYUsA9XDce9z?iOmRCb^-KYb8fE?{E`qJ5YZs5#p3V1-$3jm2#P^{n57ma!pxM?y(r%O7xb}={}QBjqN=bY)NARZ3JDGn-9(7L=Y72Az@qx;eUSVqm4 zbSgy~ps;p>w@I+6jIntq+Kk5@E}zNmb)J3s81V|IiQRpg1HEh;yra!Gu*Jrqd7@H= zfD%&Xm8KS$LD2k&6YdLQ)>B}?4~pdE6_(@k`;mfUc9(3&qn+YxrTY`-4F z>Ewuf$vL8XuTS?N@$bx6dCdOv)6*5#A5St-u3fcXQBNRTQ)f2)(+pbfS;+!hWW%wc zYbX4x+felPF%{P3!6=e&V5I!S+q2VgaSVj?8)Yk`|DZGtqRI|x*RvjaejF)T2T4#) zHleVE+;Y&x>9dRJZqIIOtlnOr$~-5{w}6pwx9dDntQ3#c{Nol1WDm1%9))*P+8|uQ z%b1*|6V6AIN+T#5lTp-`b!0QOw+4qai`Jc#_4l zf##vwHAwR?CY-S0O{&WGz6<6Hd<1mg-@;}9TPuxCh7exL8vOkCyBPNOA`UJ%G1JzG ze!&oVFEhLC*Me~>Y-4b%4D|bDP^OU<3yX%HL(hA)b8_~+N=9+evV*iCDF#a^T-rG0-@E+^(^myT~y=m_E&` zba&?4XWW%1L23E-UM5f1n}BTh4&4f^_CVP@Aq_}G(iFP41r3i~ zh+&UQm?)U|pv&*^3e=k<9Z{Aa=MRFcRLnlrX;@Uh=os@pYYP6X7x9TX`K&kqaA#X; zn|hmV*m}4-3AKMHzTp!7^C}Hgj3}|U9lzdlNIWE>4Pc=SAfvsc1{{e1j;vqIg}tl( z%VPu|if`iy+eId*_Ia4t$sEJ;mK$j`G< z2rkW2@XIe#FtAiGG_^7`w=y(RFfuSSQ6L&HzyPDCFF%ld2#71@EcJKfI%FVlXWEu7 zkB^!kQ~KkMtTyjg&G6_p)%30UP;PfN$9E#jLP7CqcO{LiXEHI&({#vtera3m@l!H> zclw0yn{1NqaqBy9=woZ{VU|-{UP^1qgdcjY_p_*cR)lo>(*pZd^Rlm8maKjBQuv1D z{l7hVqCb?<)~<|L&hk9p?%Dsy0`@}B+`kk5-`VKC`i7}m{i~P3H!kh?^ZJ1J#+9>A n`xdQTeL3;e|IPK1)2;dc2!7uxH8;7YbGf>i`BMFiGKJ3o;kbSD literal 0 HcmV?d00001 diff --git a/balloon_inputs.mat b/balloon_inputs.mat new file mode 100644 index 0000000000000000000000000000000000000000..77eda46ed0613ed74ac9a874c90501a3d7cb4fa3 GIT binary patch literal 455 zcmeZu4DoSvQZUssQ1EpO(M`+DN!3vZ$Vn_o%P-2cQV4Jk_w>_Ia4t$sEJ;mK$j`G< z2rkW2@XIe#FtAiGG_^7@w=y+XFfuSSQ6L&HzyPDCFF%m24a6059w#RxY+y(-Q+U>J zN#Ow3Q)8vX2@wG`RY`^u(h?JrCI}`ZFfbT1F!h17=)<)XplT7Y-=nb6NZsDBRru8Z zr<;#uF3+9C%uw>4eJ)6wAuert3VvraJN(52O*^&q_xv{uyPk2i0|1iQj)4FG literal 0 HcmV?d00001 diff --git a/license.txt b/license.txt new file mode 100644 index 0000000..b81e7c4 --- /dev/null +++ b/license.txt @@ -0,0 +1,24 @@ +Copyright (c) 2014, Brent Justice +All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are +met: + + * Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + * Redistributions in binary form must reproduce the above copyright + notice, this list of conditions and the following disclaimer in + the documentation and/or other materials provided with the distribution + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" +AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE +IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE +ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE +LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR +CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF +SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS +INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN +CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) +ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE +POSSIBILITY OF SUCH DAMAGE.