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 0000000..3149841 Binary files /dev/null and b/Balloon_Shape/natural_balloon.slx differ 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 0000000..e2a42b9 Binary files /dev/null and b/GUI.fig differ diff --git a/GUI.m b/GUI.m new file mode 100644 index 0000000..bbffa34 --- /dev/null +++ b/GUI.m @@ -0,0 +1,251 @@ +function varargout = GUI(varargin) +% GUI MATLAB code for GUI.fig +% GUI, by itself, creates a new GUI or raises the existing +% singleton*. +% +% H = GUI returns the handle to a new GUI or the handle to +% the existing singleton*. +% +% GUI('CALLBACK',hObject,eventData,handles,...) calls the local +% function named CALLBACK in GUI.M with the given input arguments. +% +% GUI('Property','Value',...) creates a new GUI or raises the +% existing singleton*. Starting from the left, property value pairs are +% applied to the GUI before GUI_OpeningFcn gets called. An +% unrecognized property name or invalid value makes property application +% stop. All inputs are passed to GUI_OpeningFcn via varargin. +% +% *See GUI Options on GUIDE's Tools menu. Choose "GUI allows only one +% instance to run (singleton)". +% +% See also: GUIDE, GUIDATA, GUIHANDLES + +% Edit the above text to modify the response to help GUI + +% Last Modified by GUIDE v2.5 09-Nov-2014 15:47:28 + +% Begin initialization code - DO NOT EDIT +gui_Singleton = 1; +gui_State = struct('gui_Name', mfilename, ... + 'gui_Singleton', gui_Singleton, ... + 'gui_OpeningFcn', @GUI_OpeningFcn, ... + 'gui_OutputFcn', @GUI_OutputFcn, ... + 'gui_LayoutFcn', [] , ... + 'gui_Callback', []); +if nargin && ischar(varargin{1}) + gui_State.gui_Callback = str2func(varargin{1}); +end + +if nargout + [varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:}); +else + gui_mainfcn(gui_State, varargin{:}); +end +% End initialization code - DO NOT EDIT + +% --- Executes just before GUI is made visible. +function GUI_OpeningFcn(hObject, eventdata, handles, varargin) +% Clear workspace +evalin('base','clear') +clc; + +% Dependencies +addpath(genpath(pwd)); + +% Choose default command line output for GUI +handles.output = hObject; + +% Update handles structure +guidata(hObject, handles); + +% --- Outputs from this function are returned to the command line. +function varargout = GUI_OutputFcn(hObject, eventdata, handles) +% varargout cell array for returning output args (see VARARGOUT); +% hObject handle to figure +% eventdata reserved - to be defined in a future version of MATLAB +% handles structure with handles and user data (see GUIDATA) + +% Get default command line output from handles structure +varargout{1} = handles.output; + + +% --- Executes on button press in Create_Balloon_Button. +function Create_Balloon_Button_Callback(hObject, eventdata, handles) +% hObject handle to Create_Balloon_Button (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles structure with handles and user data (see GUIDATA) + +%read input values and convert to appropriate units +if (get(handles.Hydrogen,'Value') == get(handles.Hydrogen,'Max')) + R_gas = 4124; %specific gas constant hydrogen (SI) +else + R_gas = 2077; %specific gas constant helium (SI) +end +rho_PE = str2double(get(handles.rho_PE,'String')); %kg/m^3 +thickness_PE = str2double(get(handles.thickness_PE,'String')) * 1E-6; %m +Wpayload = str2double(get(handles.Wpayload,'String')) * 4.44822162; %N +alt_apogee = str2double(get(handles.alt_apogee,'String')) * 0.3048; %m +numGores = str2double(get(handles.numGores,'String')); % +savex('balloon_inputs.mat','hObject','eventdata','handles'); %save inputs to file + +%run balloon creation code. save output to guidata +handles.balloon = create_balloon('balloon_inputs'); +guidata(hObject, handles); + + +function rho_PE_Callback(hObject, eventdata, handles) +% hObject handle to rho_PE (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles structure with handles and user data (see GUIDATA) + +% --- Executes during object creation, after setting all properties. +function rho_PE_CreateFcn(hObject, eventdata, handles) +% hObject handle to rho_PE (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles empty - handles not created until after all CreateFcns called + +% Hint: edit controls usually have a white background on Windows. +% See ISPC and COMPUTER. +if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) + set(hObject,'BackgroundColor','white'); +end + +function thickness_PE_Callback(hObject, eventdata, handles) +% hObject handle to thickness_PE (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles structure with handles and user data (see GUIDATA) + +% --- Executes during object creation, after setting all properties. +function thickness_PE_CreateFcn(hObject, eventdata, handles) +% hObject handle to thickness_PE (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles empty - handles not created until after all CreateFcns called + +% Hint: edit controls usually have a white background on Windows. +% See ISPC and COMPUTER. +if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) + set(hObject,'BackgroundColor','white'); +end + +function Wpayload_Callback(hObject, eventdata, handles) +% hObject handle to Wpayload (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles structure with handles and user data (see GUIDATA) + +% --- Executes during object creation, after setting all properties. +function Wpayload_CreateFcn(hObject, eventdata, handles) +% hObject handle to Wpayload (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles empty - handles not created until after all CreateFcns called + +% Hint: edit controls usually have a white background on Windows. +% See ISPC and COMPUTER. +if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) + set(hObject,'BackgroundColor','white'); +end + +function alt_apogee_Callback(hObject, eventdata, handles) +% hObject handle to alt_apogee (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles structure with handles and user data (see GUIDATA) + +% --- Executes during object creation, after setting all properties. +function alt_apogee_CreateFcn(hObject, eventdata, handles) +% hObject handle to alt_apogee (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles empty - handles not created until after all CreateFcns called + +% Hint: edit controls usually have a white background on Windows. +% See ISPC and COMPUTER. +if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) + set(hObject,'BackgroundColor','white'); +end + +function numGores_Callback(hObject, eventdata, handles) +% hObject handle to numGores (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles structure with handles and user data (see GUIDATA) + +% --- Executes during object creation, after setting all properties. +function numGores_CreateFcn(hObject, eventdata, handles) +% hObject handle to numGores (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles empty - handles not created until after all CreateFcns called + +% Hint: edit controls usually have a white background on Windows. +% See ISPC and COMPUTER. +if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) + set(hObject,'BackgroundColor','white'); +end + + + +% --- Executes on button press in save_gui_inputs. +function save_gui_inputs_Callback(hObject, eventdata, handles) +% hObject handle to save_gui_inputs (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles structure with handles and user data (see GUIDATA) +[filename, pathname] = uigetfile('GUI_inputs/*.mat', 'Pick or create a GUI iunputs file'); + if isequal(filename,0) || isequal(pathname,0) + % user pressed cancel. Do nothing + else + inputs = struct; + m = 1; + NAMES = fieldnames(handles); + for n = 1:1:length(NAMES) + if isprop(handles.(NAMES{n}),'style') + if strcmp('edit' , get(handles.(NAMES{n}),'style') ) + inputs(m).tag = get(handles.(NAMES{n}),'tag'); + inputs(m).string = get(handles.(NAMES{n}),'string'); + m = m+1; + end + end + end + save(fullfile(pathname,filename),'inputs'); + end + +% --- Executes on button press in load_GUI_inputs. +function load_GUI_inputs_Callback(~, eventdata, handles) +% hObject handle to load_GUI_inputs (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles structure with handles and user data (see GUIDATA) +[filename, pathname] = uigetfile('GUI_inputs/*.mat', 'Pick or create a GUI iunputs file'); + if isequal(filename,0) || isequal(pathname,0) + % user pressed cancel. Do nothing + else + temp = load(fullfile(pathname,filename)); + inputs = temp.inputs; + for n = 1:1:length(inputs) + set(handles.(inputs(n).tag),'string',inputs(n).string) + end + end + + + +function gui_inputs_filepath_Callback(hObject, eventdata, handles) +% hObject handle to gui_inputs_filepath (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles structure with handles and user data (see GUIDATA) + +% Hints: get(hObject,'String') returns contents of gui_inputs_filepath as text +% str2double(get(hObject,'String')) returns contents of gui_inputs_filepath as a double + + +% --- Executes during object creation, after setting all properties. +function gui_inputs_filepath_CreateFcn(hObject, eventdata, handles) +% hObject handle to gui_inputs_filepath (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles empty - handles not created until after all CreateFcns called + +% Hint: edit controls usually have a white background on Windows. +% See ISPC and COMPUTER. +if ispc && isequal(get(hObject,'BackgroundColor'), get(0,'defaultUicontrolBackgroundColor')) + set(hObject,'BackgroundColor','white'); +end + + +% -------------------------------------------------------------------- +function uipanel7_ButtonDownFcn(hObject, eventdata, handles) +% hObject handle to uipanel7 (see GCBO) +% eventdata reserved - to be defined in a future version of MATLAB +% handles structure with handles and user data (see GUIDATA) diff --git a/GUI_inputs/default_GUI_inputs.mat b/GUI_inputs/default_GUI_inputs.mat new file mode 100644 index 0000000..600d1f7 Binary files /dev/null and b/GUI_inputs/default_GUI_inputs.mat differ diff --git a/balloon_inputs.mat b/balloon_inputs.mat new file mode 100644 index 0000000..77eda46 Binary files /dev/null and b/balloon_inputs.mat differ 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.