Hola, estoy haciendo un simulador de la Ley de Stokes, pero no se como hacer que cuando llegue a determinada altura sepa que llego al fluido y ahí haga las ecuaciones de fluido
intente con un if, pero esto me daña el programa :S ¿Alguna sugerencia?
Gracias por la atención.
function varargout = Fluido(varargin)
% FLUIDO MATLAB code for Fluido.fig
% FLUIDO, by itself, creates a new FLUIDO or raises the existing
% singleton*.
%
% H = FLUIDO returns the handle to a new FLUIDO or the handle to
% the existing singleton*.
%
% FLUIDO('CALLBACK',hObject,eventData,handles,...) calls the local
% function named CALLBACK in FLUIDO.M with the given input arguments.
%
% FLUIDO('Property','Value',...) creates a new FLUIDO or raises the
% existing singleton*. Starting from the left, property value pairs are
% applied to the GUI before Fluido_OpeningFcn gets called. An
% unrecognized property name or invalid value makes property application
% stop. All inputs are passed to Fluido_OpeningFcn via varargin.
%
% *See GUI Options on GUIDE's Tools me nu. Choose "GUI allows only one
% instance to run (singleton)".
%
% See also: GUIDE, GUIDATA, GUIHANDLES

% Edit the above text to modify the response to help Fluido

% Last Modified by GUIDE v2.5 31-Jul-2013 16:01:53

% Begin initialization code - DO NOT EDIT
gui_Singleton = 1;
gui_State = struct('gui_Name', mfilename, ...
'gui_Singleton', gui_Singleton, ...
'gui_OpeningFcn', @Fluido_OpeningFcn, ...
'gui_OutputFcn', @Fluido_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 Fluido is made visible.
function Fluido_OpeningFcn(hObject, eventdata, handles, varargin)
% This function has no output args, see OutputFcn.
% hObject handle to figure
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)
% varargin command line arguments to Fluido (see VARARGIN)

% Choose default command line output for Fluido
handles.output = hObject;

% Update handles structure
guidata(hObject, handles);

% UIWAIT makes Fluido wait for user response (see UIRESUME)
% uiwait(handles.figure1);


% --- Outputs from this function are returned to the command line.
function varargout = Fluido_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 selection change in fluido.
function fluido_Callback(hObject, eventdata, handles)
% hObject handle to fluido (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)

% Hints: contents = cellstr(get(hObject,'String')) returns fluido contents as cell array
% contents{get(hObject,'Value')} returns selected item from fluido


val = get(hObject, 'Value');
global de;
global vis;

switch val
case 2
% El usuario seleccionó la primera opción AGUA
de=1000;
vis=0.001;
case 3
% El usuario seleccionó la segunda opción ACEITE
de=920;
vis=0.03;
case 4
% El usuario seleccionó la 3 opción GLICERINA
de=1200;
vis=1.5;
end
% --- Executes during object creation, after setting all properties.
function fluido_CreateFcn(hObject, eventdata, handles)
% hObject handle to fluido (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles empty - handles not created until after all CreateFcns called

% Hint: popupmenu 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 radio_Callback(hObject, eventdata, handles)
global r;
% hObject handle to radio (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 radio as text
% str2double(get(hObject,'String')) returns contents of radio as a double
r=(str2double(get(hObject,'String')));
% --- Executes during object creation, after setting all properties.
function radio_CreateFcn(hObject, eventdata, handles)
% hObject handle to radio (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 pushbutton1.
function pushbutton1_Callback(hObject, eventdata, handles)
% hObject handle to pushbutton1 (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles structure with handles and user data (see GUIDATA)

global r ;
r=str2double(get(handles.radio,'String'));
global de ;
global vis;

yo=15;
vyo=0;
g=-9.8;

tcaidali=(yo*2/g)^(1/2);


t=0:1/50:5:(((6)^2/7)+20);

%tMax = 0;
%deltaT = 0.50;
%t = tcaidali;



%%Ecuaciones caida libre
x=0*t;
vy = vyo + g*t;
y=yo+vyo*t+(1/2)*g*t.^2;

%%%Material--->Cobre
densidadcobre=8500;
%Entra en el fluido a una altura de cinco metros
if yo==5
%%Ecuaciones Fluidos

m=(densidadcobre*(4/3)*pi*(r.^2));
b=(6*pi*r*vis);
E=(de*(4/3)*pi*(r.^3)*g);

vl=((2*r.^2*g*(densidadcobre-de))/(9*vis));
v=vl*(1-exp((-b*t/m)));
x=vl*(t-(m/b)*(1-exp((-b/m)*t)));



for i=1:length(t) % para animar las grafica
pause(1/50) % tiempo para que se detenga tiempo igual al de arriba


%subplot(2,2,1)% Matriz 2x2 posición1
plot(handles.axes1,t(i),vy(i))
xlabel(handles.axes1,'tiempo t(s)')
ylabel(handles.axes1,'Velocidad en Y - Vy (m/s)')
axis(handles.axes1,[min(t) max(t) min(vy) max(vy)])
hold (handles.axes1,'on')
yo=yo+1;
plot(handles.axes2,x(i),y(i), 'o r')
xlabel(handles.axes2,'Posición x(m)')
ylabel(handles.axes2,'Posición y(m)')
axis(handles.axes1,axis,[-2 2 min(y) max(y)])
hold (handles.axes2,'on')

plot(handles.axes1,t(i),v(i))
xlabel(handles.axes1,'tiempo t(s)')
ylabel(handles.axes1,'Velocidad en Vl - Vl (m/s)')
axis(handles.axes1,axis,[min(t) max(t) min(v) max(v)])
hold (handles.axes1,'on')


end
end
%tMax = 3;
%deltaT = 0.01;
%t = 0;

%Bola caida libre
%while tiempoTranscurrido < tMax
%caídalibre
%x =0*deltaT;
% y = yo+Vyo*deltaT+(1/2)*g*deltaT.^2;

% vy = vy + g*deltaT;

%plot(handles.axes2, t, vy);
% xlabel(handles.axes2,'tiempo en t(s).')
% ylabel(handles.axes2,'Velocidad en Vy(m/s)')
% axis([0 tMax -5 5])
% hold (handles.axes2,'on')

%end

%v=v2*(1-exp((-b/m)*t));
%x=vt*(t-(m/b)*(1-exp((-b/m)*t)));
%vl=(2*r^2*g*(de-df))*9*vis;

%b=(6*Pi*r*vis);
%E=(df*(4/3)*Pi*(r^3)*g);
%m=(de*(4/3)*Pi*(r^2));

%plotear v vs tiempo


% --- Executes during object creation, after setting all properties.
%function axes1_CreateFcn(hObject, eventdata, handles)
% hObject handle to axes1 (see GCBO)
% eventdata reserved - to be defined in a future version of MATLAB
% handles empty - handles not created until after all CreateFcns called

% Hint: place code in OpeningFcn to populate axes1

Estuve haciendo una pequeña simulación y la fuerza de rozamiento obtenida por la expresión de Stokes sale como 4 órdenes de magnitud menos, por lo que prácticamente no afecta a la aceleración de bajada. Posiblemente el error esté en las unidades de algún dato.

Aquí mi simulación:

%Datos líquido
de=1000;
vis=1;
%Datos esfera
r=10;
densidadcobre=8500;
m=(densidadcobre*(4/3)*pi*(r.^3));
%Datos iniciales: altura inicial, velocidad inicial, aceleración y altura
%del líquido
yo=15;
vyo=0;
g=-9.8; 
yent=5;
 
%%Ecuaciones caida libre
tent=(-vyo-sqrt(vyo^2-2*g*(yo-yent)))/g;
t=0:0.001:tent;
vy = vyo + g*t;
y=yo+vyo*t+(1/2)*g*t.^2;
y1=y(length(y));
v1=vy(length(vy));
while y1>0
    b=(6*pi*r*vis*v1);
    al=(g*m+b)/m;
    v1= vyo + al*(t(length(t))+0.01);
    y1=yo+vyo*(t(length(t))+0.01)+(1/2)*al*(t(length(t))+0.01)^2;
    y=[y,y1];
    vy=[vy,v1];
    t=[t,t(length(t))+0.01];
end
 
 
 
%%Ecuaciones Fluidos
subplot(2,1,1)
plot(t,y)
xlabel('tiempo t(s)')
ylabel('Velocidad en Y - Vy (m/s)')
subplot(2,1,2)
plot(t,vy)

Eso está muy poco optimizado y es un primer ejemplo solo para el caso en el que se tira desde una altitud superior a la entrada en el líquido.