function output = pplane81(action,input1,input2,input3) % pplane81 is an interactive tool for studying planar autonomous systems of % differential equations. When pplane81 is executed, a pplane81 Setup % window is opened. The user may enter the differential % equation and specify a display window using the interactive % controls in the Setup window. Up to 4 parameters may also be % specified. In addition the user is give a choice of the type of % field displayed and the number of field points. % % When the Proceed button is pressed on the Setup window, the pplane81 % Display window is opened, and a field of the type requested is % displayed for the system. When the mouse button is depressed in % the pplane81 Display window, the solution to the system with that % initial condition is calculated and plotted. % % Other options are available in the menus. These are % fairly self explanatory. % % This is version 6.0, and will only run on version 6 of MATLAB. % Copywright (c) 1995, 1997, 1998, 1999, 2000, 2001, 2002, 2003 % John C. Polking, Rice University % Last Modified: April 24, 2003 startstr = 'pplane81'; if nargin < 1 action ='initialize'; end if get(0,'screendepth') == 1 style = 'bw'; end if strcmp(action,'initialize') % First we make sure that there is no other copy of pplane81 % running, since this causes problems. pph = findobj('tag','pplane81'); if ~isempty(pph); qstring = {'There are some pplane81 figures open although they may be invisible. ';... 'What do you want to do?'}; tstring = 'Only one copy of pplane81 can be open at one time.'; b1str = 'Restart pplane81.'; b2str = 'Just close those figures.'; b3str = 'Do nothing.'; answer = questdlg(qstring,tstring,b1str,b2str,b3str,b1str); if strcmp(answer,b1str) delete(pph); %pplane81;return eval(startstr);return elseif strcmp(answer,b2str) delete(pph);return else return end end % if ~isempty(pph); % Make sure tempdir is on the MATLABPATH. We want to be sure that we % return the path to its current position when we exit. p = path; tmpdir = tempdir; ll = length(tmpdir); tmpdir = tmpdir(1:ll-1); ud.remtd = 0; if isempty(findstr(tmpdir,p)) ud.remtd = 1; addpath(tempdir) end % Next we look for old files created by pplane81. % First in the current directory. oldfiles = dir('pptp*.m'); ofiles = cell(0,1); kk = zeros(0,1); for k = 1:length(oldfiles) fn = oldfiles(k).name; fid = fopen(fn,'r'); ll = fgetl(fid); ll = fgetl(fid); ll = fgetl(fid); fclose(fid); if strcmp(ll,'%% Created by pplane81') delete(fn) else l= length(ofiles); ofiles{l+1} = fn; end end %Then in the temp directory. oldfiles = dir([tempdir,'pptp*.m']); for k = 1:length(oldfiles) fn = [tempdir,oldfiles(k).name]; fid = fopen(fn,'r'); ll = fgetl(fid); ll = fgetl(fid); ll = fgetl(fid); fclose(fid); if strcmp(ll,'%% Created by pplane81') delete(fn) else l= length(ofiles); ofiles{l+1} = fn; end end lll = length(ofiles); if lll >0 if lll == 1 astr = 'The file'; bstr = 'If so it can be safely deleted.'; else astr = 'The files'; bstr = 'If so they can be safely deleted.'; end fprintf([astr,'\n']); for j = 1:lll fn = ofiles{j}; disp([' ',fn]); end fprintf('may have been created by an older version of PPLANE.\n'); fprintf(bstr,'\n\n'); end style = 'white'; ppdir = pwd; ssize = get(0,'defaultaxesfontsize'); npts = 20; solver = 'Dormand Prince'; tol = 1e-4; stepsize = 0.1; if exist('ppstart','file') H = ppstart; if ~isempty(H) if isfield(H,'style') style = H.style; end if isfield(H,'size') ssize = H.size; end if isfield(H,'npts') npts = H.npts; end if isfield(H,'solver') solver = H.solver; end if isfield(H,'ppdir') ppdir = H.ppdir; end if isfield(H,'stepsize') stepsize = H.stepsize; end if isfield(H,'tolerance') tol = H.tolerance; end end end if get(0,'screendepth') == 1 style = 'bw'; end ud.ssize = ssize; ud.ppdir = ppdir; comp = computer; if strcmp(comp,'PCWIN') ud.fontsize = 0.8*ud.ssize; else ud.fontsize = ud.ssize; end % Set up for the menu of systems. system.name = 'default system'; system.xvar = 'x'; system.yvar = 'y'; system.xder = ' 2*x - y + 3*(x^2-y^2) + 2*x*y'; system.yder = ' x - 3*y - 3*(x^2-y^2) + 3*x*y'; system.pname = {}; system.pval = {}; system.fieldtype = 'arrows'; system.npts = npts; system.wind = [-2 4 -4 2]; system(2).name = 'linear system'; system(2).xvar = 'x'; system(2).yvar = 'y'; system(2).xder = ' A*x + B*y'; system(2).yder = ' C*x + D*y'; system(2).pname = {'A', 'B', 'C', 'D', '', ''}; system(2).pval = {'2', '2', '-2', '-3'}; system(2).fieldtype = 'arrows'; system(2).npts = npts; system(2).wind = [-5 5 -5 5]; system(3).name = 'vibrating spring'; system(3).xvar = 'x'; system(3).yvar = 'v'; system(3).xder = ' v'; system(3).yder = ' -(k*x + d*v)/m'; system(3).pname = {'k', 'm', 'd', '', '', ''}; system(3).pval = {'3', '1', '0'}; system(3).fieldtype = 'arrows'; system(3).npts = npts; system(3).wind = [-5 5 -5 5]; system(4).name = 'pendulum'; system(4).xvar = '\theta'; system(4).yvar = '\omega'; system(4).xder = ' \omega'; system(4).yder = ' -sin(\theta) - D*\omega'; system(4).pname = {'D', '', '', '', '', ''}; system(4).pval = {'0'}; system(4).fieldtype = 'arrows'; system(4).npts = npts; system(4).wind = [-10 10 -4 4]; system(5).name = 'predator prey'; system(5).xvar = 'prey'; system(5).yvar = 'predator'; system(5).xder = ' (A - B*predator)*prey'; system(5).yder = ' (D*prey - C)*predator'; system(5).pname = {'A', 'B', 'C', 'D', '', ''}; system(5).pval = {'0.4', '0.01', '0.3', '0.005'}; system(5).fieldtype = 'arrows'; system(5).npts = npts; system(5).wind = [0 120 0 80]; system(6).name = 'competing species'; system(6).xvar = 'x'; system(6).yvar = 'y'; system(6).xder = ' r*(1 - x - A*y)*x'; system(6).yder = ' s*(1 - y - B*x)*y'; system(6).pname = {'r', 's', 'A', 'B', '', ''}; system(6).pval = {'0.4', '0.6', '5', '4'}; system(6).fieldtype = 'nullclines'; system(6).npts = npts; system(6).wind = [0 1 0 1]; system(7).name = 'cooperative species'; system(7).xvar = 'x'; system(7).yvar = 'y'; system(7).xder = ' r*(1 - x + A*y)*x'; system(7).yder = ' s*(1 - y + B*x)*y'; system(7).pname = {'r', 's', 'A', 'B', '', ''}; system(7).pval = {'0.4', '0.6', '2', '0.3'}; system(7).fieldtype = 'nullclines'; system(7).npts = npts; system(7).wind = [0 9 0 5]; system(8).name = 'van der Pol''s equation'; system(8).xvar = 'x'; system(8).yvar = 'y'; system(8).xder = ' M*x - y - x^3'; system(8).yder = ' x'; system(8).pname = {'M', '', '', '', '', ''}; system(8).pval = {'2'}; system(8).fieldtype = 'arrows'; system(8).npts = npts; system(8).wind = [-5 5 -5 5]; system(9).name = 'Duffing''s equation'; system(9).xvar = 'x'; system(9).yvar = 'y'; system(9).xder = ' y'; system(9).yder = ' -(k*x + c*y + l*x^3)/m'; system(9).pname = {'k', 'c', 'm', 'l', '', ''}; system(9).pval = {'-1', '0', '1', '1'}; system(9).fieldtype = 'arrows'; system(9).npts = npts; system(9).wind = [-3 3 -3 3]; system(10).name = 'square limit set'; system(10).xvar = 'x'; system(10).yvar = 'y'; system(10).xder = ' (y + x/5)*(1-x^2)'; system(10).yder = ' -x*(1-y^2)'; system(10).pname = {'', '', '', '', '', ''}; system(10).pval = {'', '', '', ''}; system(10).fieldtype = 'arrows'; system(10).npts = npts; system(10).wind = [-1.5 1.5 -1 1]; ud.c = system(1); % Changed values. pname = ud.c.pname; for kk = length(pname)+1:6 pname{kk} = ''; end pval = ud.c.pval; for kk = length(pval)+1:6 pval{kk} = ''; end ud.c.pname = pname; ud.c.pval = pval; ud.o = ud.c; % Original values. % ud.h = system(1); % This will be the handles in the % setup window. ud.style = style; ud.npts = npts; ud.settings.magn = 1.25; ud.settings.refine = 8; ud.settings.tol = tol; ud.settings.solver = solver; ud.settings.stepsize = stepsize; ud.settings.speed = 100; ud.system = system; ud.solvers = {'ppdp45'; 'rk4'; 'ode45'; 'ode23'; 'ode113'; 'ode15s'; 'ode23s'; 'ode23t'; 'ode23tb'}; switch style case 'black' color.temp = [1 0 0]; % red for temporary orbits color.orb = [1 1 0]; % yellow for orbits color.eqpt = [1 0 0]; % red for eq. pts. % color.arrows = [0 1 1]; % cyan for arrows % color.arrows = [.5 .5 .9]; % purple for arrows color.arrows = .5*[1 1 1]; % gray for arrows color.narrows = .7*[1 1 1]; % gray for nullcline arrows color.tx = [1 1 0]; % yellow for xt plots & 3D plots color.ty = [1 0 0]; % red for yt plots color.ta = [1 0 0]; % red for axis plots color.sep = [.2,1,0]; % green for separatrices color.xcline = [1 1 1]; % white for xclines color.ycline = [1 0 1]; % magenta for yclines color.level = [1,.5,.5]; case 'white' color.temp = [1 0 0]; % red for temporary orbits color.orb = [0 0 1]; % blue for orbits color.eqpt = [1 0 0]; % red for eq. pts. color.arrows = 0.7*[1 1 1]; % gray for arrows color.narrows = .4*[1 1 1]; % gray for nullcline arrows color.tx = [0 0 1]; % blue for xt plots & 3D plots color.ty = [1 0 0]; % red for yt plots color.ta = [1 0 0]; % red for axis plots color.sep = [0,1,0];% green for separatrices color.xcline = [1 0 .75]; % magenta for xclines color.ycline = [1 .5 0]; % orange for yclines color.level = 0.8*[.9,.5,.8]; case 'test' color.temp = [1 0 0]; % red for temporary orbits color.orb = [0 0 1]; % blue for orbits color.arrows = .7*[1 1 1]; % gray for arrows color.eqpt = [1 0 0]; % red for eq. pts. color.narrows = .4*[1 1 1]; % gray for nullcline arrows color.tx = [0 0 1]; % blue for xt plots & 3D plots color.ty = [1 0 0]; % red for yt plots color.ta = [1 0 0]; % red for axis plots color.sep = [0,1,0];% green for separatrices color.xcline = [1 0 .75]; % magenta for xclines color.ycline = [1 .5 0]; % orange for yclines color.level = [1,.5,.5]; case 'display' color.temp = [1 0 0]; % red for temporary orbits color.orb = [0 0 1]; % blue for orbits color.arrows = .4*[1 1 1]; % gray for arrows color.eqpt = [1 0 0]; % red for eq. pts. color.narrows = .4*[1 1 1]; % gray for nullcline arrows color.tx = [0 0 1]; % blue for xt plots & 3D plots color.ty = [1 0 0]; % red for yt plots color.ta = [1 0 0]; % red for axis plots color.sep = 2*[.5 0 .5];% [0,1,0];% green for separatrices color.xcline = [1 0 .75]; % magenta for xclines color.ycline = [1 .5 0]; % orange for yclines color.level = 0.8*[.9,.5,.8]; case 'bw' color.temp = [1 1 1]; % white for everything color.orb = [1 1 1]; color.eqpt = [1 1 1]; color.arrows = [1 1 1]; color.narrows = [1 1 1]; color.tx = [1 1 1]; color.ty = [1 1 1]; color.ta = [1 1 1]; color.sep = [1 1 1]; color.xcline = [1 1 1]; color.ycline = [1 1 1]; color.level = [1,1,1]; end ud.color = color; ud.level = ' '; ppset = figure('name','pplane81 Setup','NumberTitle','off',... 'tag','pplane81','visible','off',... 'UserData',ud); pplane81('figdefault',ppset); frame(1) = uicontrol('style','frame','visible','off'); eq(1)=uicontrol('style','text',... 'horizon','center',... 'string','The differential equations.','visible','off'); ext = get(eq(1),'extent'); rr=ext(4)/10; texth =ext(4)+4; % 19; % Height of text boxes. varw = 40*rr; % Length of variable boxes. equalw =13*rr; % Length of equals.(30) eqlength = 230*rr; % Length of right hand sides of equations. winstrlen = 120*rr; % Length of string boxes in display frame. left = 0; % Left margin of the frames. frsep = 1; %3; % Separation between frames. separation = texth; % Separation between boxes. dfigwidth =2*left + varw+equalw+eqlength+10; % Width of the figure. dfigurebot = 30; % Bottom of the figure. buttw = dfigwidth/3; qwind = [0,frsep,buttw,separation]; % Quit button rwind = [buttw,frsep,buttw,separation]; % Revert " pwind = [2*buttw,frsep,buttw,separation]; % Proceed " disfrbot = 2*frsep + separation; % Display frame. disfrw = winstrlen + varw +10; disfrht = 5*separation + 10; disfrwind = [left, disfrbot, disfrw, disfrht]; pfrbot = disfrbot + disfrht +frsep; % Parameter frame. pfrw = dfigwidth -2*left; pfrht = 3*separation + 10; pfrwind = [left, pfrbot, pfrw, pfrht]; defrbot = pfrbot + pfrht + frsep; % Equation frame. defrw = pfrw; defrht = 3*separation + 10; defrwind = [left, defrbot, defrw, defrht]; ffrbot = disfrbot; % Field frame. ffrleft = left + disfrw + frsep; ffrw = dfigwidth -left - ffrleft; ffrht = disfrht; ffrwind = [ffrleft, ffrbot, ffrw, ffrht]; dfigureheight = defrbot + defrht +frsep; % Height of the figure. set(ppset,'pos',[30 dfigurebot dfigwidth dfigureheight]); set(frame(1),'pos',defrwind); xname=[ 'ud = get(gcf,''UserData'');'... 'Xname=get(ud.h.xvar,''string'');'... 'minxstr = [''The minimum value of '',Xname,'' = ''];',... 'set(ud.h.twind(1),''string'',minxstr);'... 'maxxstr = [''The maximum value of '',Xname,'' = ''];',... 'set(ud.h.twind(2),''string'',maxxstr);'... 'ud.c.xvar = Xname;'... 'ud.flag = 0;'... 'ud.c.name = '''';',... 'set(gcf,''UserData'',ud);']; yname=[ 'ud = get(gcf,''UserData'');'... 'Yname=get(ud.h.yvar,''string'');'... 'minystr = [''The minimum value of '',Yname,'' = ''];',... 'set(ud.h.twind(3),''string'',minystr);'... 'maxystr = [''The maximum value of '',Yname,'' = ''];',... 'set(ud.h.twind(4),''string'',maxystr);'... 'ud.c.yvar = Yname;'... 'ud.flag = 0;'... 'ud.c.name = '''';',... 'set(gcf,''UserData'',ud);']; xder =[ 'ud = get(gcf,''UserData'');'... 'ud.c.xder = get(ud.h.xder,''string'');'... 'ud.flag = 0;'... 'ud.c.name = '''';',... 'set(gcf,''UserData'',ud);']; yder =[ 'ud = get(gcf,''UserData'');'... 'ud.c.yder = get(ud.h.yder,''string'');'... 'ud.flag = 0;'... 'ud.c.name = '''';',... 'set(gcf,''UserData'',ud);']; equationbot = defrbot + 5; eqlabelbot = equationbot + 2*separation; xbot = equationbot + separation; % Bottom of x equation. ybot = equationbot; % Bottom of y equation. lablen =200*rr; eqlableft = (dfigwidth-lablen)/2; eqleft = left + 5; fudge = 0.15*separation; set(eq(1),'pos',[eqlableft eqlabelbot lablen texth]); tcolor = get(gcf,'defaultuicontrolbackgroundcolor'); ecolor = 'w'; ud.h.xvar=uicontrol('pos',[eqleft, xbot, varw, texth],... 'style','edit',... 'horizon','right',... 'string',ud.o.xvar,... 'call',xname,... 'backgroundcolor',ecolor,... 'visible','off'); eq(2) = uicontrol('style','text',... 'pos',[eqleft+varw xbot equalw texth],... 'horizon','center',... 'string',''' = ',... 'backgroundcolor',tcolor,... 'visible','off'); ud.h.xder=uicontrol('pos',[eqleft+varw+equalw xbot eqlength texth],... 'string',ud.o.xder,... 'horizon','left','style','edit',... 'backgroundcolor',ecolor,... 'call',xder,'visible','off'); ud.h.yvar=uicontrol('pos',[eqleft ybot varw texth],... 'style','edit',... 'horizon','right',... 'string',ud.o.yvar,... 'backgroundcolor',ecolor,... 'call',yname,'visible','off'); eq(3) = uicontrol('style','text',... 'pos',[eqleft+varw ybot equalw texth],... 'horizon','center','string',''' = ',... 'backgroundcolor',tcolor',... 'visible','off'); ud.h.yder=uicontrol('pos',[eqleft+varw + equalw ybot eqlength texth],... 'string',ud.o.yder,... 'horizon','left','style','edit',... 'backgroundcolor',ecolor,... 'call',yder,'visible','off'); frame(2) = uicontrol('style','frame','pos',disfrwind,'visible','off'); w1 = [ 'ud = get(gcf,''UserData'');'... 'nnn = str2num(get(ud.h.wind(1),''string''));'... 'if isempty(nnn),',... ' set(ud.h.wind(1),''string'',''?'');',... ' nnn = NaN;',... 'end,',... 'ud.c.wind(1) = nnn;',... 'ud.c.name = '''';',... 'set(gcf,''UserData'',ud);']; w2 = [ 'ud = get(gcf,''UserData'');'... 'nnn = str2num(get(ud.h.wind(2),''string''));'... 'if isempty(nnn),',... ' set(ud.h.wind(2),''string'',''?'');',... ' nnn = NaN;',... 'end,',... 'ud.c.wind(2) = nnn;',... 'ud.c.name = '''';',... 'set(gcf,''UserData'',ud);']; w3 = [ 'ud = get(gcf,''UserData'');'... 'nnn = str2num(get(ud.h.wind(3),''string''));'... 'if isempty(nnn),',... ' set(ud.h.wind(3),''string'',''?'');',... ' nnn = NaN;',... 'end,',... 'ud.c.wind(3) = nnn;',... 'ud.c.name = '''';',... 'set(gcf,''UserData'',ud);']; w4 = [ 'ud = get(gcf,''UserData'');'... 'nnn = str2num(get(ud.h.wind(4),''string''));'... 'if isempty(nnn),',... ' set(ud.h.wind(4),''string'',''?'');',... ' nnn = NaN;',... 'end,',... 'ud.c.wind(4) = nnn;',... 'ud.c.name = '''';',... 'set(gcf,''UserData'',ud);']; winbot1 = disfrbot + disfrht - 5 - separation; winbot2 = winbot1 - separation; winbot3 = winbot2 - separation; winbot4 = winbot3 - separation; winbot5 = winbot4 - separation; dwindow = uicontrol('style','text',... 'pos',[eqleft winbot1 disfrw-10 texth],... 'horizon','center',... 'string','The display window.','visible','off'); % ud.h.twind contains the handles to the text windows, and ud.h.wind % contains the handles to the edit windows. twstr1 = ['The minimum value of ',ud.o.xvar,' = ']; ud.h.twind(1) = uicontrol('style','text',... 'pos',[eqleft winbot2 winstrlen texth],... 'horizon','right',... 'string',twstr1,... 'backgroundcolor',tcolor',... 'visible','off'); ud.h.wind(1) = uicontrol('style','edit',... 'pos',[eqleft+winstrlen winbot2 40*rr texth],... 'string',num2str(ud.o.wind(1)),... 'backgroundcolor',ecolor,... 'call',w1,'visible','off'); twstr2 = ['The maximum value of ',ud.o.xvar,' = ']; ud.h.twind(2) = uicontrol('style','text',... 'pos',[eqleft winbot3 winstrlen texth],... 'horizon','right',... 'string',twstr2,... 'backgroundcolor',tcolor',... 'visible','off'); ud.h.wind(2) = uicontrol('style','edit',... 'pos',[eqleft+winstrlen winbot3 40*rr texth],... 'string',num2str(ud.o.wind(2)),... 'backgroundcolor',ecolor,... 'call',w2,'visible','off'); twstr3 = ['The minimum value of ',ud.o.yvar,' = ']; ud.h.twind(3)= uicontrol('style','text',... 'pos',[eqleft winbot4 winstrlen texth],... 'horizon','right',... 'string',twstr3,... 'backgroundcolor',tcolor',... 'visible','off'); ud.h.wind(3) = uicontrol('style','edit',... 'pos',[eqleft+winstrlen winbot4 40*rr texth],... 'string',num2str(ud.o.wind(3)),... 'backgroundcolor',ecolor,... 'call',w3,'visible','off'); twstr4 = ['The maximum value of ',ud.o.yvar,' = ']; ud.h.twind(4) = uicontrol('style','text',... 'pos',[eqleft winbot5 winstrlen texth],... 'horizon','right',... 'string',twstr4,... 'backgroundcolor',tcolor',... 'visible','off'); ud.h.wind(4) = uicontrol('style','edit',... 'pos',[eqleft+winstrlen winbot5 40*rr texth],... 'string',num2str(ud.o.wind(4)),... 'backgroundcolor',ecolor,... 'call',w4,'visible','off'); frame(3)=uicontrol('style','frame','pos',pfrwind,'visible','off'); pncall = [ '[h,fig] = gcbo;'... 'ud = get(fig,''UserData'');'... 'num = get(h,''UserData'');'... 'ud.c.pname{num} = get(ud.h.pname(num),''string'');'... 'ud.flag = 0;'... 'set(gcf,''UserData'',ud);']; pvcall = [ '[h,fig] = gcbo;'... 'ud = get(fig,''UserData'');'... 'num = get(h,''UserData'');'... 'ud.c.pval{num} = get(ud.h.pval(num),''string'');'... 'ud.flag = 0;'... 'set(gcf,''UserData'',ud);']; pnamew = 40*rr; pvalw = 50*rr; peqw = 10*rr; pbot(3) = pfrbot + 5; pbot(2) = pbot(3) + separation; pbot(1) = pbot(2) + separation; pleft1 = eqleft + 50*rr + 5; peqleft1 = pleft1 + pnamew; pvleft1 = peqleft1 + peqw; pleft2 = dfigwidth - 10 - pnamew - pvalw - peqw; peqleft2 = pleft2 + pnamew; pvleft2 = peqleft2 + peqw; paratit=uicontrol('style','text',... 'horizon','center',... 'string',{'Parameters';'or';'expressions'},... 'backgroundcolor',tcolor',... 'visible','off'); ext = get(paratit,'extent'); paratitw = ext(3); pos = [eqleft pfrbot+2+texth/2 paratitw 2.1*texth]; set(paratit,'pos',pos); psep = 20; pvalw = (dfigwidth - 2*eqleft - paratitw)/2 - psep - pnamew - peqw; pval = ud.c.pval; pname = ud.c.pname; for jj = 1:3 for kk = 1:2 pleft = eqleft + paratitw +psep +(kk-1)*(pnamew+peqw+pvalw+ psep); peqleft = pleft + pnamew; pvleft = peqleft + peqw; K = kk +2*(jj-1); name = pname{K}; value = pval{K}; ud.h.pname(K) = uicontrol('style','edit',... 'pos',[pleft pbot(jj) pnamew texth],... 'horizon','right','string',name,... 'UserData',K,... 'call',pncall,... 'visible','off',... 'backgroundcolor','w'); equal(K) = uicontrol('style','text',... 'pos',[peqleft pbot(jj)-fudge peqw texth],... 'horizon','center',... 'string','=',... 'visible','off'); ud.h.pval(K) = uicontrol('style','edit',... 'pos',[pvleft pbot(jj) pvalw texth],... 'string',value,... 'call',pvcall,... 'visible','off',... 'UserData',K,... 'backgroundcolor','w'); end end ud.c.pname = pname; ud.c.pval = pval; butt(1) = uicontrol('style','push',... 'pos',qwind,... 'string','Quit','call',... 'pplane81(''quit'')',... 'visible','off'); butt(2) = uicontrol('style','push',... 'pos',rwind,... 'string','Revert',... 'call','pplane81(''revert'')',... 'visible','off'); butt(3) = uicontrol('style','push',... 'pos',pwind,... 'string','Proceed',... 'call','pplane81(''proceed'')',... 'visible','off'); fframe = uicontrol('style','frame','pos',ffrwind,'visible','off'); ffrtitle = uicontrol('style','text',... 'pos',[ffrleft+5,winbot1,ffrw-10,texth],... 'string','The direction field.',... 'horizon','center','visible','off'); radleft = ffrleft + 3; radw = 50*rr; typewindw = radw +6; typewind = [ffrleft, ffrbot, typewindw, ffrht-separation-3]; textwindl = ffrleft+typewindw; textleft = textwindl + 3; textw = ffrw - typewindw; textwind = [textwindl, ffrbot,textw, ffrht-separation-3]; typeframe = uicontrol('style','frame','pos',typewind,'visible','off'); textframe = uicontrol('style','frame','pos',textwind,'visible','off'); switch ud.o.fieldtype case 'nullclines' rval1 = 1;rval2 = 0;rval3 = 0;rval4 = 0; case 'lines' rval1 = 0;rval2 = 2;rval3 = 0;rval4 = 0; case 'arrows' rval1 = 0;rval2 = 0;rval3 = 3;rval4 = 0; case 'none' rval1 = 0;rval2 = 0;rval3 = 0;rval4 = 4; otherwise error(['Unknown fieldtype ',ud.o.fieldtype,'.']) end ud.h.rad(1) = uicontrol('style','radio',... 'pos',[radleft winbot4 radw texth],... 'string','Nullclines',... 'value',rval1,... 'visible','off'); ud.h.rad(2) = uicontrol('style','radio',... 'pos',[radleft winbot3 radw texth],... 'string','Lines',... 'value',rval2,... 'max',2,... 'visible','off'); ud.h.rad(3) = uicontrol('style','radio',... 'pos',[radleft winbot2 radw texth],... 'string','Arrows',... 'value',rval3,... 'max',3,... 'visible','off'); ud.h.rad(4) = uicontrol('style','radio',... 'pos',[radleft winbot5 radw texth],... 'string','None',... 'value',rval4,... 'max',4,... 'visible','off'); for i=1:4 set(ud.h.rad(i),'UserData',ud.h.rad(:,[1:(i-1),(i+1):4])); end callrad = [ 'me = get(gcf,''currentobject'');',... 'kk = get(me,''max'');',... 'set(get(me,''UserData''),''value'',0),',... 'set(me,''value'',kk);',... 'ud = get(gcf,''UserData'');',... 'switch kk,',... ' case 1, ud.c.fieldtype = ''nullclines'';',... ' case 2, ud.c.fieldtype = ''lines'';',... ' case 3, ud.c.fieldtype = ''arrows'';',... ' case 4, ud.c.fieldtype = ''none'';',... 'end,',... 'set(gcf,''UserData'',ud);']; set(ud.h.rad,'call',callrad); nfptsstr = {'Number of'; 'field points per'; 'row or column.'}; nfptstext = uicontrol('style','text',... 'pos',[textleft winbot4 textw-5 2.5*texth],... 'string',nfptsstr,... 'horizon','center',... 'visible','off'); callnfpts = [ 'ppset = findobj(''name'',''pplane81 Setup'');',... 'ud = get(ppset,''UserData'');'... 'me = ud.h.npts;',... 'kk = str2num(get(me,''string''));',... 'if isempty(kk),',... ' set(me,''string'',''?'');',... ' kk = NaN;',... 'else,',... ' kk = floor(kk);',... ' [m,N] = computer;'... ' if (N <= 8192),',... ' N = 32;',... ' else,',... ' N = 50;',... ' end,'... ' kk = min([N,max([5,kk])]);'... ' set(me,''string'',num2str(kk));'... 'end,'... 'ud.c.npts = kk;',... 'set(ppset,''UserData'',ud);']; npos = [textleft+(textw -30*rr)/2,winbot5 30*rr,texth]; ud.h.npts = uicontrol('style','edit',... 'pos',npos,... 'string',ud.o.npts,... 'call',callnfpts,... 'backgroundcolor','w',... 'visible','off'); delgall = ['sud = get(gcf,''UserData'');',... 'mh = get(sud.h.gallery,''children'');',... 'add = findobj(sud.h.gallery,''tag'',''add system'');',... 'mh(find(mh == add)) = [];',... 'delete(mh);',... 'set(sud.h.gallery,''UserData'',[]);',... 'set(findobj(''tag'',''load default''),''enable'',''on'')']; % Menus hhsetup = get(0,'showhiddenhandles'); set(0,'showhiddenhandles','on'); mefile = findobj(ppset,'label','&File'); meedit = findobj(ppset,'label','&Edit'); delete(findobj(ppset,'label','&Tools')); delete(findobj(ppset,'label','&View')); delete(findobj(ppset,'label','&Insert')); % File menu meexp = findobj(mefile,'label','&Export...'); meprev = findobj(mefile,'label','Print Pre&view...'); mepset = findobj(mefile,'label','Pa&ge Setup...'); set(get(mefile,'child'),'vis','off'); meload = uimenu(mefile,'label','Load a system ...',... 'call','pplane81(''loadsyst'',''system'');',... 'pos',1); mesave = uimenu(mefile,'label','Save the current system ...',... 'call','pplane81(''savesyst'',''system'');',... 'pos',2); meloadg = uimenu(mefile,'label','Load a gallery ...',... 'call','pplane81(''loadsyst'',''gallery'');',... 'separator','on','pos',3); mesaveg = uimenu(mefile,'label','Save a gallery ...',... 'call','pplane81(''savesyst'',''gallery'');',... 'tag','savegal','pos',4); medelg = uimenu(mefile,'label','Delete the current gallery',... 'call',delgall,'pos',5); melddg = uimenu(mefile,'label','Load the default gallery',... 'call','pplane81(''loadsyst'',''default'');',... 'enable','on',... 'tag','load default','pos',6); merevert = uimenu(mefile,'label','Revert','call',... 'pplane81(''revert'')',... 'separator','on','pos',7); meproceed = uimenu(mefile,... 'label','Proceed',... 'call','pplane81(''proceed'')',... 'separator','off',... 'accelerator','G','pos',8); set(mepset,'vis','on','pos',9); set(meprev,'vis','on','pos',10); set(meexp,'vis','on','pos',11,'separator','off'); merestart = uimenu(mefile,'label',... 'Restart pplane81',... 'call','pplane81(''restart'')',... 'separator','on','pos',12); mequit = uimenu(mefile,... 'label','Quit pplane81',... 'call','pplane81(''quit'')',... 'separator','off','pos',13); % Edit menu set(get(meedit,'child'),'vis','off'); meclrf = uimenu(meedit,'label','Clear equations',... 'call',['ud = get(gcf,''UserData'');h = ud.h;',... 'set([h.xvar,h.xder,h.yvar,h.yder],''string'','''');'],... 'accelerator','E'); pclear = [ 'ud = get(gcf,''UserData'');h = ud.h;',... 'set([h.pname,h.pval],''string'','''');',... 'ud.c.pname = {'''','''','''','''','''','''',''''};',... 'ud.c.pval = {'''','''','''','''','''','''',''''};',... 'set(gcf,''UserData'',ud);',... ]; meclrp = uimenu(meedit,'label','Clear parameters',... 'call',pclear,... 'accelerator','N'); meclrwind = uimenu(meedit,'label','Clear display window',... 'call',['ud = get(gcf,''UserData'');',... 'set(ud.h.wind,''string'','''');'],... 'accelerator','D'); allclear = [ 'ud = get(gcf,''UserData'');h = ud.h;',... 'set([h.xvar,h.xder,h.yvar,h.yder],''string'','''');',... 'set([h.pname,h.pval,h.wind],''string'','''');',... 'ud.c.pname = {'''','''','''','''','''','''',''''};',... 'ud.c.pval = {'''','''','''','''','''','''',''''};',... 'set(gcf,''UserData'',ud);',... ]; meclrall = uimenu(meedit,'label','Clear all',... 'call',allclear,... 'accelerator','A',... 'separator','on'); % Gallery menu. sysmenu = uimenu('label','Gallery','visible','off','pos',3); meadd = uimenu(sysmenu,'label','Add current system to the gallery',... 'call','pplane81(''addgall'');','tag','add system'); sep = 'on'; for kk = 1:length(system) kkk = num2str(kk); if kk == 2, sep = 'off';end sysmen(kk) = uimenu(sysmenu,'label',system(kk).name,... 'call',['pplane81(''system'',',kkk,')'],... 'separator',sep,'visible','off'); end set(sysmenu,'UserData',system); ud.h.gallery = sysmenu; ud.flag = 0; ud.egg = (exist('EASTEREGG') ==2); % Record the handles in the User Data of the Set Up figure. set(ppset,'UserData',ud); hhhh = findobj(ppset,'type','uicontrol'); set(hhhh,'units','normal') set(ppset,'visible','on','resize','on'); set(get(ppset,'children'),'visible','on'); set(get(sysmenu,'children'),'visible','on'); %end set(0,'showhiddenhandles',hhsetup); elseif strcmp(action,'savesyst') ppset = findobj('name','pplane81 Setup'); type = input1; sud = get(ppset,'UserData'); switch type case 'system' systems = get(sud.h.gallery,'UserData'); newsyst = sud.c; fn = newsyst.name; if ~isempty(fn) fn(find(abs(fn)==32))='_'; % Replace spaces by underlines. end fn = [fn, '.pps']; comp = computer; switch comp case 'PCWIN' filter = [sud.ppdir, '\',fn]; case 'MAC2' filter = [sud.ppdir,':', fn]; otherwise filter = [sud.ppdir,'/', fn]; end [fname,pname] = uiputfile(filter,'Save the system as:'); if fname == 0,return;end if ~strcmp(fname,fn) ll = length(fname); if (ll>4 & strcmp(fname(ll-3:ll),'.pps')) fn = fname; else fn = [fname, '.pps']; end newsyst.name = fn; sud.c.name = fn; set(ppset,'UserData',sud); end newsysts = newsyst; case 'gallery' systems = get(sud.h.gallery,'UserData'); ll = length(systems); if ll == 0 warndlg(['There are no systems to make up a gallery.'],'Warning'); return end names = cell(ll,1); for j=1:ll names{j} = systems(j).name; end [sel,ok] = listdlg('PromptString','Select the systems',... 'Name','Gallery selection',... 'ListString',names); if isempty(sel) return else newsysts = systems(sel); end comp = computer; switch comp case 'PCWIN' prompt = [sud.ppdir,'\*.ppg']; case 'MAC2' prompt = [sud.ppdir,':*.ppg']; otherwise prompt = [sud.ppdir,':/.ppg']; end [fname,pname] = uiputfile(prompt,'Save the gallery as:'); ll = length(fname); if (ll>4 & strcmp(fname(ll-3:ll),'.ppg')) fn = fname; else fn = [fname, '.ppg']; end newsyst.name = fn(1:ll-4); sud.c.name = fn(1:ll-4); set(ppset,'UserData',sud); end % switch type ll = length(newsysts); fid = fopen([pname fn],'w'); ppstring = '%%%% PPLANE file %%%%'; fprintf(fid,[ppstring,'\n']); for k = 1:ll fprintf(fid,'\n'); nstr = newsysts(k).name; nstr = strrep(nstr,'''',''''''); nstr = ['H.name = ''', nstr, ''';\n']; fprintf(fid,nstr); xname = newsysts(k).xvar; xnstr = ['H.xvar = ''', xname]; xnstr = strrep(xnstr,'\','\\'); xnstr = [xnstr, ''';\n']; fprintf(fid,xnstr); yname = newsysts(k).yvar; ynstr = ['H.yvar = ''', yname]; ynstr = strrep(ynstr,'\','\\'); ynstr = [ynstr, ''';\n']; fprintf(fid,ynstr); xder = newsysts(k).xder; xdstr = ['H.xder = ''', xder]; xdstr = strrep(xdstr,'\','\\'); xdstr = [xdstr, ''';\n']; fprintf(fid,xdstr); yder = newsysts(k).yder; ydstr = ['H.yder = ''', yder]; ydstr = strrep(ydstr,'\','\\'); ydstr = [ydstr, ''';\n']; fprintf(fid,ydstr); pname = strrep(newsysts(k).pname,'\','\\'); pval = strrep(newsysts(k).pval,'\','\\'); pnl = length(pname); pvl = length(pval); for kk = 1:6 if kk <= pnl pns = pname{kk}; else pns = ''; end if kk <= pvl pvs = pval{kk}; else pvs = ''; end if kk == 1 pnstr = ['H.pname = {''', pns, '''']; pvstr = ['H.pval = {''', pvs, '''']; else pnstr = [pnstr, ',''',pns, '''']; pvstr = [pvstr, ',''',pvs, '''']; end end pnstr = [pnstr, '};\n']; pvstr = [pvstr, '};\n']; fprintf(fid,pnstr); fprintf(fid,pvstr); ftstr = ['H.fieldtype = ''', newsysts(k).fieldtype]; ftstr - strrep(ftstr,'\','\\'); ftstr = [ftstr, ''';\n']; fprintf(fid,ftstr); nstr = ['H.npts = ', num2str(newsysts(k).npts),';\n']; fprintf(fid,nstr); wind = newsysts(k).wind; wstr = ['H.wind = [', num2str(wind),'];\n']; fprintf(fid,wstr); end fclose(fid); elseif strcmp(action,'loadsyst') % This loads either a system or a gallery. sud = get(gcf,'UserData'); pos = get(gcf,'pos'); wpos = [pos(1),pos(2)+pos(4)+20,300,20]; waith = figure('pos',wpos,... 'NumberTitle','off',... 'vis','off',... 'next','replace',... 'menubar','none',... 'resize','off',... 'createfcn',''); axes('pos',[0.01,0.01,0.98,0.98],... 'vis','off'); xp = [0 0 0 0]; yp = [0 0 1 1]; xl = [1 0 0 1 1]; yl = [0 0 1 1 0]; patchh = patch(xp,yp,'r','edgecolor','r','erase','none'); lineh = line(xl,yl,'erase','none','color','k'); type = input1; set(sud.h.gallery,'enable','off'); if strcmp(type,'default') set(waith,'name','Loading the default gallery.','vis','on'); set(findobj('tag','load default'),'enable','off'); megall = sud.h.gallery; mh = get(megall,'children'); add = findobj(megall,'tag','add system'); mh(find(mh == add)) = []; delete(mh); newsysstruct = get(megall,'UserData'); system = sud.system; ll = length(system); x = 1/(ll+2); xp = [xp(2),x,x,xp(2)]; set(patchh,'xdata',xp); set(lineh,'xdata',xl); drawnow; sep = 'on'; for kk = 1:length(system) kkk = num2str(kk); if kk ==2, sep = 'off';end uimenu(megall,'label',system(kk).name,... 'call',['pplane81(''system'',',kkk,')'],... 'separator',sep); end % for set(megall,'UserData',system); else comp = computer; switch comp case 'PCWIN' prompt = [sud.ppdir,'\']; case 'MAC2' prompt = [sud.ppdir,':']; otherwise prompt = [sud.ppdir,'/']; end if strcmp(type,'system') prompt = [prompt,'*.pps']; [fname,pname] = uigetfile(prompt,'Select a system to load.'); elseif strcmp(type,'gallery') prompt = [prompt,'*.ppg']; [fname,pname] = uigetfile(prompt,'Select a gallery to load.'); end % if strcmp if fname == 0 delete(waith); set(sud.h.gallery,'enable','on'); return; end set(waith,'name',['Loading ',fname],'vis','on'); fid = fopen([pname fname],'r'); sline = fgetl(fid); if strcmp(sline,'%% PPLANE file %%') date = 'new'; else date = 'old'; end newsysts = {}; switch date case 'old' newsysts{1} = sline; kk = 1; case 'new' kk = 0; end while ~feof(fid) kk = kk + 1; newsysts{kk} = fgetl(fid); end fclose(fid); newsysts = newsysts([1:kk]); false = 0; switch date case 'old' if mod(kk,19) false = 1; end case 'new' if mod(kk,11) false = 1; end end %switch date if false if strcmp(type,'system') warndlg(['The file ',fname, ' does not define a proper system.'],... 'Warning'); elseif strcmp(type,'gallery') warndlg(['The file ',fname, ' does not define a proper gallery.'],... 'Warning'); end set(sud.h.gallery,'enable','on'); delete(waith); return end %if false switch date case 'old' x = 19/(kk+38); xp = [xp(2),x,x,xp(2)]; set(patchh,'xdata',xp); set(lineh,'xdata',xl); drawnow; nnn = kk/19; flds = fieldnames(sud.c); flds=flds(:); for j = 0:(nnn-1) newsystemp = newsysts([(j*19+1):(j+1)*19]); newsyst.name = newsystemp{1}; newsyst.xvar = newsystemp{2}; newsyst.yvar = newsystemp{3}; newsyst.xder = newsystemp{4}; newsyst.yder = newsystemp{5}; newsyst.pname = {newsystemp{6}, newsystemp{7},... newsystemp{8}, newsystemp{9},'',''}; newsyst.pval = {newsystemp{10},... newsystemp{11},... newsystemp{12},... newsystemp{13},'',''}; newsyst.fieldtype = newsystemp{14}; newsyst.npts = str2num(newsystemp{15}); wind = newsystemp(16:19); newsyst.wind = [str2num(wind{1}),str2num(wind{2}),... str2num(wind{3}),str2num(wind{4})]; newsysstruct(j+1) = newsyst; end % for j = 0:(nnn-1) case 'new' x = 11/(kk+22); xp = [xp(2),x,x,xp(2)]; set(patchh,'xdata',xp); set(lineh,'xdata',xl); drawnow; nnn = kk/11; for j = 1:nnn for k = 2:11; eval(newsysts{(j-1)*11+k}); end newsysstruct(j) = H; end end %switch date end % if strcmp(type,'default') & else nnn = length(newsysstruct); ignoresyst = {}; for j = 1:nnn x = (j+1)/(nnn+2); xp = [xp(2),x,x,xp(2)]; set(patchh,'xdata',xp); set(lineh,'xdata',xl); drawnow; newsyst = newsysstruct(j); sname = newsyst.name; sname(find(abs(sname) == 95)) = ' '; % Replace underscores with spaces. newsyst.name = sname; ignore = pplane81('addgall',newsyst); if ignore == -1; ignoresyst{length(ignoresyst)+1} = sname; end end % for j = 1:nnn l = length(ignoresyst); if l % There was at least one system which was a dup with a % different name. if l == 1 message = {['The system ',ignoresyst{1},'" duplicates a ',... 'system already in the gallery and was not added.']}; else message = 'The systems '; for k = 1:(l-1) message = [message,'"',ignoresyst{k},'", ']; end message = {[message,'and "',ignoresyst{l},'" duplicate ',... 'systems already in the gallery and were not added.']}; end % if l == 1 & else helpdlg(message,'Ignored systems'); end % if l if strcmp(type,'system') % Added a system. if ignore > 0 % The system was ignored. kk = ignore; else systems = get(sud.h.gallery,'UserData'); kk = length(systems); end pplane81('system',kk); end if strcmp('type','default') pplane81('system',1); end set(sud.h.gallery,'enable','on'); x = 1; xp = [xp(2),x,x,xp(2)]; set(patchh,'xdata',xp); set(lineh,'xdata',xl); drawnow; delete(waith); elseif strcmp(action,'addgall') output = 0; ppset = findobj('name','pplane81 Setup'); sud = get(ppset,'UserData'); if nargin < 2 % We are adding the current system. syst = sud.c; snstr = 'Provide a name for this system.'; sname = inputdlg(snstr,'System name',1,{syst.name}); if isempty(sname),return;end sname = sname{1}; if ~strcmp(sname,syst.name) sud.c.name = sname; set(ppset,'UserData',sud); syst.name = sname; end else % We have a system coming from a file. syst = input1; sname = syst.name; end pnl = length(syst.pname); for kk = (pnl+1):6 syst.pname{kk} = ''; end pvl = length(syst.pval); for kk = (pvl+1):6 syst.pval{kk} = ''; end systems = get(sud.h.gallery,'UserData'); ll = length(systems); kk = 1; while ((kk<=ll) & (~strcmp(sname,systems(kk).name))) kk = kk + 1; end nameflag = (kk<=ll); ssyst = rmfield(syst,'name'); kk = 1; while ((kk<=ll) & (~isequal(ssyst,rmfield(systems(kk),'name')))) kk = kk + 1; end systflag = 2*(kk<=ll); flag = nameflag + systflag; switch flag case 1 % Same name but different system. mh = findobj(sud.h.gallery,'label',sname); prompt = {['The system "',sname,'", which you wish to ',... 'add to the gallery has ',... 'the same name as a different system ',... 'already in the gallery. Please ',... 'specify the name for the newly added system.'],... 'Specify the name for the old system.'}; title = 'Two systems with the same name'; lineno = 1; defans = {sname,sname}; answer = inputdlg(prompt,title,lineno,defans); if isempty(answer),return,end sname = answer{1}; systems(kk).name = answer{2}; set(mh,'label',answer{2}); output = kk; case 2 % Two names for the same system. oldname = systems(kk).name; mh = findobj(sud.h.gallery,'label',oldname); prompt = {['The system "',sname,'", which you wish to add ',... 'to the gallery is the same as a system which is ',... 'already in the gallery with the name "',oldname,'". ',... 'Please specify which name you wish to use.']}; title = 'Two names for the same system.'; lineno = 1; defans = {oldname}; answer = inputdlg(prompt,title,lineno,defans); if isempty(answer),return,end systems(kk).name = answer{1}; set(mh,'label',answer{1}); output = kk; case 3 % Systems and names the same. output = -1; otherwise end % switch set(sud.h.gallery,'UserData',systems); syst.name = sname; if flag <=1 switch ll case 0 systems = syst; sepstr = 'on'; case 10 systems(11) = syst; if strcmp(systems(10).name,'square limit set') sepstr = 'on'; else sepstr = 'off'; end otherwise systems(ll+1) = syst; sepstr = 'off'; end kkk = num2str(ll+1); newmenu = uimenu(sud.h.gallery,'label',sname,... 'call',['pplane81(''system'',',kkk,')'],... 'separator',sepstr); set(findobj('tag','savegal'),'enable','on'); end set(sud.h.gallery,'UserData',systems); elseif strcmp(action,'system') ppset = findobj('name','pplane81 Setup'); ud = get(ppset,'UserData'); kk = input1; if isstr(kk) kk = str2num(input1); end system = get(ud.h.gallery,'UserData'); syst = system(kk); xname = syst.xvar; yname = syst.yvar; set(ud.h.xvar,'string',xname); set(ud.h.yvar,'string',yname); set(ud.h.xder,'string',syst.xder); set(ud.h.yder,'string',syst.yder); pname = syst.pname; pval = syst.pval; pnl = length(pname); pvl = length(pval); for kk = 1:6 if kk <= pnl set(ud.h.pname(kk),'string',pname{kk}); else set(ud.h.pname(kk),'string',''); syst.pname{kk} = ''; end if kk <= pvl set(ud.h.pval(kk),'string',pval{kk}); else set(ud.h.pval(kk),'string',''); syst.pval{kk} = ''; end end ud.o = syst; ud.c = syst; set(ud.h.twind(1),'string',['The minimum value of ',xname,' = ']); set(ud.h.twind(2),'string',['The maximum value of ',xname,' = ']); set(ud.h.twind(3),'string',['The minimum value of ',yname,' = ']); set(ud.h.twind(4),'string',['The maximum value of ',yname,' = ']); for kk = 1:4 set(ud.h.wind(kk),'string',num2str(syst.wind(kk))); end set(ud.h.npts,'string',num2str(syst.npts)); switch syst.fieldtype case 'nullclines' rval = [1 0 0 0]; case 'lines' rval = [0 2 0 0]; case 'arrows' rval = [0 0 3 0]; case 'none' rval = [0 0 0 4]; otherwise error(['Unknown fieldtype ',ud.o.fieldtype,'.']) end for i=1:4 set(ud.h.rad(i),'value',rval(i)); end ud.flag = 0; set(ppset,'UserData',ud); elseif strcmp(action,'revert') ud = get(gcf,'UserData'); ud.c = ud.o; syst = ud.o; xname = syst.xvar; yname = syst.yvar; set(ud.h.xvar,'string',xname); set(ud.h.yvar,'string',yname); set(ud.h.xder,'string',syst.xder); set(ud.h.yder,'string',syst.yder); pname = syst.pname; pval = syst.pval; pnl = length(pname); pvl = length(pval); for kk = 1:6 if kk <= pnl nstr = pname(kk); else nstr = ''; end if kk <= pvl vstr = pval(kk); else vstr = ''; end set(ud.h.pname(kk),'string',nstr); set(ud.h.pval(kk),'string',vstr); end set(ud.h.twind(1),'string',['The minimum value of ',xname,' = ']); set(ud.h.twind(2),'string',['The maximum value of ',xname,' = ']); set(ud.h.twind(3),'string',['The minimum value of ',yname,' = ']); set(ud.h.twind(4),'string',['The maximum value of ',yname,' = ']); for kk = 1:4 set(ud.h.wind(kk),'string',num2str(syst.wind(kk))); end set(ud.h.npts,'string',num2str(syst.npts)); switch syst.fieldtype case 'lines' rval(1) = 1; rval(2) = 0;rval(3) = 0; case 'arrows' rval(1) = 0;rval(2) = 2;rval(3) = 0; case 'none' rval(1) = 0;rval(2) = 0;rval(3) = 3; otherwise error(['Unknown fieldtype ',ud.o.fieldtype,'.']) end for i=1:3 set(ud.h.rad(i),'value',rval(i)); end set(gcf,'UserData',ud); elseif strcmp(action,'proceed') % Proceed connects Setup with the Display window. ppset = gcf; sud = get(ppset,'UserData'); sud.o = sud.c; set(ppset,'UserData',sud); he = findobj('name','pplane81 Equilibrium point data'); hl = findobj('name','pplane81 Linearization'); close([he;hl]); % set([he;hl],'vis','off'); % Some error checking that has to be done no matter what. WINvect = sud.c.wind; if any(isnan(WINvect)) sud.flag = 0; set(ppset,'UserData',sud); errmsg = ['One of the entries defining the display window ',... 'is not a number.']; fprintf('\a') errordlg(errmsg,'PPLANE error','on'); return end xstr = sud.c.xvar; if isempty(xstr) sud.flag = 0; set(ppset,'UserData',sud); errmsg = 'The first dependent variable needs a name.'; fprintf('\a') errordlg(errmsg,'PPLANE error','on'); return end ystr = sud.c.yvar; if isempty(ystr) sud.flag = 0; set(ppset,'UserData',sud); errmsg = 'The second dependent variable needs a name.'; fprintf('\a') errordlg(errmsg,'PPLANE error','on'); return end if WINvect(2)<= WINvect(1) sud.flag = 0; set(ppset,'UserData',sud); errmsg = ['The minimum value of ', xstr,... ' must be smaller than the maximum value.']; fprintf('\a') errordlg(errmsg,'PPLANE error','on'); return end if WINvect(4)<= WINvect(3) sud.flag = 0; set(ppset,'UserData',sud); errmsg = ['The minimum value of ', ystr,... ' must be smaller than the maximum value.']; fprintf('\a') errordlg(errmsg,'PPLANE error','on'); return end if isnan(sud.c.npts) sud.flag = 0; set(ppset,'UserData',sud); errmsg = 'The entry for the number of field points is not a number.'; fprintf('\a') errordlg(errmsg,'PPLANE error','on'); return end % sud.flag = 0 if this is the first time through for this equation, % sud.flag = 1 if only the window dimensions or the field data % have been changed. % If sud.flag == 1 we only have to update things. if (sud.flag == 1) Arrflag = sud.c.fieldtype; NumbFPts = sud.c.npts; ppdisp = findobj('name','pplane81 Display'); dud = get(ppdisp,'UserData'); aud = get(dud.axes,'UserData'); wind = sud.c.wind(:); if (~all(wind == dud.syst.wind(:))) dwind = [wind(1); wind(3); -wind(2); -wind(4)]; aud.DY = [wind(2)-wind(1); wind(4)-wind(3)]; aud.cwind = dwind - dud.settings.magn*[aud.DY;aud.DY]; set(dud.axes,'UserData',aud); end arr = dud.arr; menull = findobj('tag','null'); switch Arrflag case 'nullclines' % set([arr.hx;arr.hy;arr.barrows],'vis','on'); set([arr.hx;arr.hy],'vis','on'); set([arr.lines;arr.arrows],'vis','off'); set(menull,'enable','on','label','Hide nullclines.'); case 'lines' set(arr.lines,'vis','on'); set([arr.hx;arr.hy;arr.arrows;arr.barrows],'vis','off'); set(menull,'enable','on','label','Show nullclines.'); case 'arrows' set(arr.arrows,'vis','on'); set([arr.lines;arr.hx;arr.hy;arr.barrows],'vis','off'); set(menull,'enable','on','label','Show nullclines.'); otherwise set([arr.hx;arr.hy;arr.lines;arr.arrows;arr.barrows],'vis','off'); set(menull,'enable','on','label','Show nullclines.'); end dud.syst.fieldtype = Arrflag; set(ppdisp,'UserData',dud); if ( (NumbFPts ~= dud.syst.npts) | (any(WINvect ~= dud.syst.wind) ) ) dud.syst.wind = WINvect; dud.syst.npts = NumbFPts; set(ppdisp,'UserData',dud); pplane81('dirfield',ppdisp); end figure(ppdisp); else sud.flag = 1; set(ppset,'UserData',sud); sud = get(ppset,'UserData'); % WINvect = sud.c.wind; Xname = sud.c.xvar; Yname = sud.c.yvar; xderivstr = sud.c.xder; yderivstr = sud.c.yder; pname = sud.c.pname; parav = sud.c.pval; % Convert the parameters to their current values. First remove the % blanks. Also remove the periods inserted by users attempting to % make the function array smart. xderivstr(find(abs(xderivstr)==32))=[]; l=length(xderivstr); for ( k = fliplr(findstr('.',xderivstr))) if (find('*/^' == xderivstr(k+1))) xderivstr = [xderivstr(1:k-1), xderivstr(k+1:l)] end l=l-1; end yderivstr(find(abs(yderivstr)==32))=[]; l=length(yderivstr); for ( k = fliplr(findstr('.',yderivstr))) if (find('*/^' == yderivstr(k+1))) yderivstr = [yderivstr(1:k-1), yderivstr(k+1:l)]; end l=l-1; end for kk = 1:6 pval = parav{kk}; if ~isempty(pval) pabs = abs(pval); kkk = find(pabs==32); pval(kkk) = []; l = length(pval); for ( k = fliplr(findstr('.',pval))) if (find('*/^' == pval(k+1))) pval = [pval(1:k-1), pval(k+1:l)]; end l=l-1; end parav{kk} = pval; end end % Build strings for the title. txderstr = xderivstr; tyderstr = yderivstr; kxder = find(abs(txderstr)==42); txderstr(kxder)=' '*ones(size(kxder)); txderstr = strrep(txderstr,'-',' - '); txderstr = strrep(txderstr,'+',' + '); if (abs(txderstr(1)) == 32) txderstr = txderstr(2:length(txderstr)); end kxder = find(abs(tyderstr)==42); tyderstr(kxder)=' '*ones(size(kxder)); tyderstr = strrep(tyderstr,'-',' - '); tyderstr = strrep(tyderstr,'+',' + '); if (abs(tyderstr(1)) == 32) tyderstr = tyderstr(2:length(tyderstr)); end tstr1 = [Xname,' '' = ', txderstr]; tstr2 = [Yname,' '' = ', tyderstr]; tstr = str2mat(tstr1,tstr2); dud.tstr = tstr; pstr1 = {' ';' '}; pstr2 = {' ';' '}; pstr3 = {' ';' '}; pstring = cell(6,1); for kk = 1:6 if ~isempty(parav{kk}) tpstr = parav{kk}; kxder = find(abs(tpstr)==42); % Get rid of *s tpstr(kxder)=' '*ones(size(kxder)); tpstr = strrep(tpstr,'-',' - '); % Extra spaces tpstr = strrep(tpstr,'+',' + '); if (abs(tpstr(1)) == 32) % Get rid of starting space tpstr = tpstr(2:length(tpstr)); end pstring{kk} = [pname{kk},' = ', tpstr]; else % pstring{kk} = [pname{kk},' = ',parav{kk}]; pstring{kk} = [pname{kk},' = ']; end end % Get ready to do some error trapping. SS = warning; warning off XxXxXx = WINvect(1) + rand*(WINvect(2)-WINvect(1)); YyYyYy = WINvect(3) + rand*(WINvect(4)-WINvect(3)); err = 0; % Now we remove the backslashes (\) used to get Greek into the % labels. txname = Xname; tyname = Yname; xderivstr(find(abs(xderivstr)==92))=[]; yderivstr(find(abs(yderivstr)==92))=[]; Xname(find(abs(Xname)==92))=[]; Yname(find(abs(Yname)==92))=[]; eval([Xname,'=XxXxXx;'],'err = 1;'); if err sud.flag = 0; set(ppset,'UserData',sud); errmsg = ['"',xstr, '" is not a valid variable name in MATLAB.']; fprintf('\a') errordlg(errmsg,'PPLANE error','on'); return end err = 0; eval([Yname,'=YyYyYy;'],'err = 1;'); if err sud.flag = 0; set(ppset,'UserData',sud); errmsg = ['"',ystr, '" is not a valid variable name in MATLAB.']; fprintf('\a') errordlg(errmsg,'PPLANE error','on'); return end % Replace the parameters/expressions with their values. pflag = zeros(1,6); perr = []; for kk = 1:6 if ~isempty(pname{kk}) pn = pname{kk}; pv = parav{kk}; if isempty(pv) perr = [perr, sud.h.pval(kk)]; else if isempty(str2num(pv)) % This is an expression. tpv = pv; tpv(find(abs(tpv)==92))=[]; err = 0; pval = ''; eval(['pval = ',tpv,';'],'err=1;'); if (err | isempty(pval)) errmsg = ['The expression for ',pn,' is not valid.']; fprintf('\a') errordlg(errmsg,'pplane81 error','on'); warning(SS) return end end xxderivstr = pplane81('paraeval',pn,pv,xderivstr); yyderivstr = pplane81('paraeval',pn,pv,yderivstr); if (~strcmp(xxderivstr,xderivstr) | ~strcmp(yyderivstr,yderivstr) ) pflag(kk) = 1; xderivstr = xxderivstr; yderivstr = yyderivstr; end end end end % We have to make the derivative strings array smart. l = length(xderivstr); for (k=fliplr(find((xderivstr=='^')|(xderivstr=='*')|(xderivstr=='/')))) xderivstr = [xderivstr(1:k-1) '.' xderivstr(k:l)]; l = l+1; end l = length(yderivstr); for (k=fliplr(find((yderivstr=='^')|(yderivstr=='*')|(yderivstr=='/')))) yderivstr = [yderivstr(1:k-1) '.' yderivstr(k:l)]; l = l+1; end % Some more error trapping. err = 0;res = 1; eval(['res = ',xderivstr, ';'],'err = 1;'); if err | isempty(res) if isempty(perr) errmsg = ['The first differential equation ',... 'is not entered correctly.']; else errstr1 = ['The first differential equation ',... 'does not evaluate correctly.']; errstr2 = ['At least one of the parameter values is not ',... 'a number.']; errmsg = str2mat(errstr1,errstr2); perr set(perr,'string','?'); end sud.flag = 0; set(ppset,'UserData',sud); fprintf('\a') errordlg(errmsg,'PPLANE error','on'); return; end err = 0;res = 1; eval(['res = ',yderivstr, ';'],'err = 1;'); if err | isempty(res) if isempty(perr) errmsg = ['The second differential equation ',... 'is not entered correctly.']; else errstr1 = ['The second differential equation ',... 'does not evaluate correctly.']; errstr2 = ['At least one of the parameter values is not ',... 'a number.']; errmsg = str2mat(errstr1,errstr2); set(perr,'string','?'); end sud.flag = 0; set(ppset,'UserData',sud); fprintf('\a') errordlg(errmsg,'PPLANE error','on'); return; end % If an old function m-file exists delete it, and then build a new one. % if (~strcmp(dfcn,'') & exist(dfcn)==2) delete([dfcn,'.m']);end tee = clock; tee = ceil(tee(6)*100); dfcn=['pptp',num2str(tee)]; fcnstr = ['function YYyYypr = ',dfcn,'(t,YyYy)\n\n']; commstr = '%%%% Created by pplane81\n\n'; varstr = [Xname,' = YyYy(1,:);', Yname,' = YyYy(2,:);\n\n']; lenstr = ['l = length(YyYy(1,:));\n']; derstr1 = ['XxXxxpr = ', xderivstr,';\n']; derstr2 = ['if (length(XxXxxpr) < l) XxXxxpr = XxXxxpr*ones(1,l);end\n']; derstr3 = ['YyYyypr = ', yderivstr,';\n']; derstr4 = ['if (length(YyYyypr) < l) YyYyypr = YyYyypr*ones(1,l);end\n']; derstr5 = 'YYyYypr = [XxXxxpr;YyYyypr];\n'; ppf = fopen([tempdir,dfcn,'.m'],'w'); fprintf(ppf,fcnstr); fprintf(ppf,commstr); fprintf(ppf,varstr); fprintf(ppf,lenstr); fprintf(ppf,derstr1); fprintf(ppf,derstr2); fprintf(ppf,derstr3); fprintf(ppf,derstr4); fprintf(ppf,derstr5); fclose(ppf); % Find pplane81 Display if it exists. % If pplane81 Display exists, update it. If it does not build it. ppdisp = findobj('name','pplane81 Display'); if (~isempty(ppdisp)) figure(ppdisp); dud = get(ppdisp,'UserData'); dud.syst = sud.c; dud.settings = sud.settings; dfcnn = dud.function; if exist(dfcnn)==2 delete([tempdir,dfcnn,'.m']); end xmstr = [txname,' vs. t']; ymstr = [tyname,' vs. t']; set(dud.menu(3),'label',xmstr); set(dud.menu(5),'label',ymstr); menull = findobj('tag','null'); if ~isempty(menull) delete(get(menull,'UserData')); set(menull,'UserData',[]); end else ppdisp = figure('name','pplane81 Display',... 'NumberTitle','off',... 'interrupt','on',... 'visible','off',... 'tag','pplane81'); pplane81('figdefault',ppdisp); dud = get(ppdisp,'UserData'); dud.syst = sud.c; switch dud.syst.name case 'pendulum' dud.level = 'omega^2 - 2*cos(theta)'; otherwise dud.level = ' '; end dud.settings = sud.settings; dud.egg = sud.egg; dud.noticeflag = 1; dud.contours = zeros(0,1); fs = dud.fontsize; ssize = dud.ssize; r = ssize/10; ppaxw = 437*1.2; % Default axes width ppaxh = 315*1.2; % Default axes height ppaxl = 45*1.2; % Default axes left buttw = 40*1.2; % Default button width titleh = 45; % Default title height. This is changed later. nframeh = 70; % Default notice frame height ppaxb = 4+nframeh+35; bottomedge = 38; ppdh = bottomedge + nframeh + ppaxh + titleh; uni = get(0,'units'); set(0,'units','pixels'); ss = get(0,'screensize'); set(0,'units',uni); sw = ss(3);sh = ss(4); bottom = 10; if r*ppdh > sh -bottom -35; r = (sh-bottom-35)/ppdh; fs = 10*r; lw = 0.5*r; set(gcf,'defaultaxesfontsize',fs,'defaultaxeslinewidth',lw); set(gcf,'defaulttextfontsize',fs); set(gcf,'defaultlinelinewidth',lw); set(gcf,'defaultuicontrolfontsize',fs*0.9); end % Set up the bulletin window. nframe = uicontrol('style','frame','visible','on'); nstr = {'More';'than';'five';'lines';'of text'}; dud.notice = uicontrol('style','text',... 'horiz','left',... 'string',nstr,'visible','on'); ext = get(dud.notice,'extent'); nframeh = ext(4)+2; titleh = r*titleh; ppaxl = r*ppaxl; ppaxw = r*ppaxw; ppaxh = r*ppaxh; ppaxb = nframeh+r*bottomedge; buttw = r*buttw; butth = fs+10*r; buttl = ppaxl + ppaxw + 5; buttsep = (ppaxh - butth)/2; % Set up the coordinate display cstr = '(0.99999999, 0,99999999)'; dud.ccwind = uicontrol('style','text',... 'horiz','left',... 'string',cstr,... 'visible','on'); cext = get(dud.ccwind,'extent'); ccwindtxt = uicontrol('style','text',... 'horiz','left',... 'string','Cursor position: ',... 'visible','on'); cwh = cext(4); cww = cext(3); % Set up the plot axes. dud.axes = axes('units','pix',... 'position',[ppaxl,ppaxb,ppaxw,ppaxh],... 'next','add',... 'box','on',... 'interrupt','on',... 'xgrid','on',... 'ygrid','on',... 'drawmode','fast',... 'visible','off',... 'tag','display axes'); % Set up the title. dud.title.axes = axes('box','off','xlim',[0 1],'ylim',[0 1],... 'units','pix','vis','off',... 'xtick',[-1],'ytick',[-1],... 'xticklabel','','yticklabel',''); dud.title.eq = text(0.01,0.5,' ','vert','middle'); dud.title.p1 = text(0.75,0.5,' ','vert','middle'); dud.title.p2 = text(0.65,0.5,' ','vert','middle'); dud.title.p3 = text(0.55,0.5,' ','vert','middle'); tstr = {'x_2';'y^2'}; set(dud.title.eq,'string',tstr,'units','pix'); ext = get(dud.title.eq ,'extent'); titleh = ext(4)+15*r; set(dud.title.eq,'units','data'); taxpos = [ppaxl,ppaxb+ppaxh,ppaxw,titleh]; set(dud.title.axes,'pos',taxpos,'color',get(gcf,'color')); % Finish the positions. ppdw = buttl + buttw +5; ppdh = ppaxb+ppaxh+titleh; set(nframe,'pos',[10,1,ppdw-20,nframeh]); set(dud.notice,'pos',[15,1,ppdw-30,nframeh-2],... 'string',{' ';' ';' ';' ';' '}); ctext = get(ccwindtxt,'extent'); cc1pos = [ppaxl,2+nframeh,ctext(3),cwh]; cc2pos = [ppaxl+ctext(3),2+nframeh,cww,cwh]; set(ccwindtxt,'pos',cc1pos); set(dud.ccwind,'pos',cc2pos,... 'string',''); ppdleft = max((sw-ppdw)/2,sw-ppdw-60); ppdbot = sh - ppdh - 35; ppdpos = [ppdleft,ppdbot,ppdw,ppdh]; set(ppdisp,'resize','on'); set(ppdisp,'pos',ppdpos); Arrflag = sud.c.fieldtype; % Set up the buttons stopstr = 'aud = get(gca,''UserData'');aud.stop = 4;set(gca,''UserData'',aud);'; dbutt(1) = uicontrol('style','push',... 'pos',[buttl,ppaxb+2*buttsep,buttw,butth],... 'string','Stop','call',stopstr,... 'vis','off','tag','stop'); dbutt(2) = uicontrol('style','push',... 'pos',[buttl,ppaxb,buttw,butth],... 'string','Quit',... 'call','pplane81(''quit'')','visible','off'); dbutt(3) = uicontrol('style','push',... 'pos',[buttl,ppaxb+buttsep,buttw,butth],... 'string','Print',... 'call','pplane81(''print'')','visible','off'); dud.butt = dbutt; % Menus and Toolbar hhsetup = get(0,'showhiddenhandles'); set(0,'showhiddenhandles','on'); % Configure the Toolbar. fixtb = ['set(gcbo,''state'',''off'');']; set(ppdisp,'ToolBar','none'); % Menus tmenu = findobj(ppdisp,'label','&Tools'); delete(tmenu); % Insert menu imenu = findobj(gcf,'label','&Insert'); inschild = get(imenu,'child'); legitem = findobj(inschild,'label','&Legend'); colitem = findobj(inschild,'label','&Colorbar'); delete([legitem,colitem]); % File menu fmenu = findobj(ppdisp,'label','&File',... 'pos',1); delete(findobj(fmenu,'label','&New Figure')); delete(findobj(fmenu,'label','&Open...')); delete(findobj(fmenu,'label','&Close')); set(findobj(fmenu,'label','&Save'),... 'pos',1,'separator','off'); set(findobj(fmenu,'label','Save &As...'),... 'pos',2); set(findobj(fmenu,'label','&Export...'),... 'pos',3); delete(findobj(fmenu,'label','Pre&ferences...')); set(findobj(fmenu,'label','Pa&ge Setup...'),'pos',4); set(findobj(fmenu,'label','Print Set&up...'),'pos',5); set(findobj(fmenu,'label','Print Pre&view...'),'pos',6); set(findobj(fmenu,'label','&Print...'),'pos',7); merestart = uimenu(fmenu,'label',... 'Restart pplane81',... 'call','pplane81(''restart'')',... 'separator','on'); mequit = uimenu(fmenu,'label','Quit pplane81',... 'call','pplane81(''quit'')','separator','off'); % Edit menu emenu = findobj(ppdisp,'label','&Edit',... 'pos',2); menu(2) = uimenu(emenu,'label','Zoom in.',... 'call','pplane81(''zoomin'')',... 'pos',1); zsqmenu = uimenu(emenu,'label','Zoom in square.',... 'call','pplane81(''zoominsq'')',... 'pos',2); zbmenu = uimenu(emenu,'label','Zoom back.','call',... 'pplane81(''zoomback'')',... 'enable','off','tag','zbmenu',... 'pos',3); medallsol = uimenu(emenu,'label','Erase all solutions.',... 'call','pplane81(''dallsol'')',... 'separator','on','pos',4); medallep = uimenu(emenu,'label','Erase all equilibrium points.',... 'call','pplane81(''dallep'')',... 'separator','off',... 'pos',5); medallics = uimenu(emenu,'label','Erase all marked initial points.',... 'call','pplane81(''dallics'')',... 'separator','off',... 'pos',6); medalllev = uimenu(emenu,'label','Erase all level curves.',... 'call','pplane81(''dalllev'')',... 'separator','off',... 'pos',7); medall = uimenu(emenu,'label','Erase all graphics objects.',... 'call','pplane81(''dall'')',... 'separator','off',... 'pos',8); medel = uimenu(emenu,'label','Delete a graphics object.',... 'call','pplane81(''delete'')',... 'visible','on',... 'pos',9); menutext = uimenu(emenu,... 'label','Enter text on the Display Window.',... 'call','pplane81(''text'')',... 'separator','on',... 'pos',10); set(findobj(emenu,'label','&Undo'),'separator','on',... 'pos',11); set(findobj(emenu,'label','Cu&t'),'pos',12); set(findobj(emenu,'label','&Copy'),'pos',13); set(findobj(emenu,'label','&Paste'),'pos',14); set(findobj(emenu,'label','Clea&r'),'pos',15); set(findobj(emenu,'label','&Select All'),'pos',16); set(findobj(emenu,'label','Copy &Figure'),'pos',17); set(findobj(emenu,'label','Copy &Options'),'pos',18); set(findobj(emenu,'label','F&igure Properties'),'pos',19); set(findobj(emenu,'label','&Axes Properties'),'pos',20); set(findobj(emenu,'label','C&urrent Object Properties'),'pos',21); % Graph menu megraph = uimenu('label','Graph','visible','off','pos',4); menu(3) = uimenu(megraph,'label',[txname,' vs. t'],... 'call','pplane81(''plotxy'',1)'); menu(5) = uimenu(megraph,'label',[tyname,' vs. t'],... 'call','pplane81(''plotxy'',2)'); meplot3 = uimenu(megraph,'label','Both',... 'call','pplane81(''plotxy'',3)'); meplot4 = uimenu(megraph,'label','3 D',... 'call','pplane81(''plotxy'',4)'); meplot5 = uimenu(megraph,'label','Composite',... 'call','pplane81(''plotxy'',5)'); % Solutions menu solmenu = uimenu('label','Solutions','pos',3); menukey = uimenu(solmenu,'label','Keyboard input.','call',... 'pplane81(''kbd'')'); mesev = uimenu(solmenu,'label','Plot several solutions.',... 'call','pplane81(''several'')'); meeqpt = uimenu(solmenu,'label','Find an equilibrium point.',... 'call','pplane81(''eqpt'')','separator','on'); menu(4) = uimenu(solmenu,... 'label','List computed equilibrium points.',... 'call','pplane81(''eqptlist'')'); mestunst= uimenu(solmenu,... 'label','Plot stable and unstable orbits.',... 'call','pplane81(''stunst'')','interrupt','on'); meperiod = uimenu(solmenu,... 'label', 'Find a nearly closed orbit'); periodstr = ['ud = get(gcf,''UserData'');',... 'me = gcbo;',... 'ud.dir = get(me,''UserData'');',... 'set(gcf,''UserData'',ud);',... 'pplane81(''periodic'');']; dud.period(1) = uimenu(meperiod,... 'label','forward',... 'UserData',1,... 'call', periodstr); dud.period(2) = uimenu(meperiod,... 'label','backward',... 'UserData',-1,... 'call', periodstr); dud.period(3) = uimenu(meperiod,... 'label','in both directions',... 'UserData',0,... 'call', periodstr); nullcall = ['me = gcbo;',... 'dud = get(gcf,''UserData'');',... 'handx = [dud.arr.hx; dud.arr.hy];',... 'handr = [dud.arr.hr; dud.arr.hth];',... 'arron = get(dud.arr.arrows,''vis'');',... 'lab = get(me,''label'');',... 'switch lab,',... ' case ''Show nullclines.'',',... ' set(handx,''vis'',''on'');',... ' if strcmp(arron,''on''),',... ' set(dud.arr.barrows,''vis'',''off'');',... ' else,',... ' set(dud.arr.barrows,''vis'',''on'');',... ' end,',... ' set(handr,''vis'',''off'');',... ' set(me,''label'',''Hide nullclines.'');',... ' rnh = findobj(''tag'',''rnull'');'... ' set(rnh,''label'',''Show polar nullclines.'');',... ' case ''Hide nullclines.'','... ' set(handx,''vis'',''off'');',... ' set(dud.arr.barrows,''vis'',''off'');',... ' set(me,''label'',''Show nullclines.'');',... 'end']; rnullcall = ['me = gcbo;',... 'dud = get(gcf,''UserData'');',... 'handx = [dud.arr.hx; dud.arr.hy];',... 'handr = [dud.arr.hr; dud.arr.hth];',... 'lab = get(me,''label'');',... 'switch lab,',... ' case ''Show polar nullclines.'',',... ' set(handr,''vis'',''on'');',... ' set(handx,''vis'',''off'');',... ' set(dud.arr.barrows,''vis'',''off'');',... ' set(me,''label'',''Hide polar nullclines.'');',... ' nh = findobj(''tag'',''null'');'... ' set(nh,''label'',''Show nullclines.'');',... ' case ''Hide polar nullclines.'','... ' set(handr,''vis'',''off'');',... ' set(me,''label'',''Show polar nullclines.'');',... 'end']; menunull = uimenu(solmenu,'label','Show nullclines.',... 'call',nullcall,'separator','on','tag','null'); menulevel = uimenu(solmenu,'label','Plot level curves.',... 'call','pplane81(''level'')',... 'separator','off','tag','level'); if dud.egg menurnull = uimenu(solmenu,'label','Show polar nullclines.',... 'call',rnullcall,'separator','off','tag','rnull'); metest = uimenu(solmenu,'label','Test case',... 'call','pplane81(''test case'')',... 'separator','on','tag','testcase'); end % Options menu menu(1) = uimenu('label','Options','visible','off','pos',5); meset = uimenu(menu(1),'label','Settings.',... 'call','pplane81(''settings'')'); mesolve = uimenu(menu(1),'label','Solver.'); solverstr = ['ud = get(gcf,''UserData'');',... 'me = gcbo;',... 'meud = get(me,''UserData'');',... 'ud.settings.refine = meud.refine;',... 'ud.settings.tol = meud.tol;',... 'ud.settings.solver = meud.solver;',... 'ud.settings.stepsize = meud.stepsize;',... 'set(ud.solver,''checked'',''off'');',... 'set(me,''checked'',''on'');',... 'set(gcf,''UserData'',ud);',... 'ppset = findobj(''name'',''pplane81 Setup'');',... 'sud = get(ppset,''UserData'');',... 'sud.settings = ud.settings;',... 'set(ppset,''UserData'',sud);',... 'pplane81(''settings'');']; solver = dud.settings.solver; dpset.refine = 8; dpset.tol = dud.settings.tol; dpset.solver = 'Dormand Prince'; dpset.stepsize = dud.settings.stepsize; dpset.hmax = 0; if strcmp(solver,'Dormand Prince') dpch = 'on'; else dpch = 'off'; end rk4set.refine = 1; rk4set.tol = dud.settings.tol; rk4set.solver = 'Runge-Kutta 4'; rk4set.stepsize = dud.settings.stepsize; rk4set.hmax = 0; if strcmp(solver,'Runge-Kutta 4') rk4ch = 'on'; else rk4ch = 'off'; end ode15sset.refine = 1; ode15sset.tol = dud.settings.tol; ode15sset.solver = 'ode15s'; ode15sset.stepsize = dud.settings.stepsize; ode15sset.hmax = 0; if strcmp(solver,'ode15s') ode15sch = 'on'; else ode15sch = 'off'; end ode23sset.refine = 1; ode23sset.tol = dud.settings.tol; ode23sset.solver = 'ode23s'; ode23sset.stepsize = dud.settings.stepsize; ode23sset.hmax = 0; if strcmp(solver,'ode23s') ode23sch = 'on'; else ode23sch = 'off'; end ode113set.refine = 1; ode113set.tol = dud.settings.tol; ode113set.solver = 'ode113'; ode113set.stepsize = dud.settings.stepsize; ode113set.hmax = 0; if strcmp(solver,'ode113') ode113ch = 'on'; else ode113ch = 'off'; end ode23set.refine = 1; ode23set.tol = dud.settings.tol; ode23set.solver = 'ode23'; ode23set.stepsize = dud.settings.stepsize; ode23set.hmax = 0; if strcmp(solver,'ode23') ode23ch = 'on'; else ode23ch = 'off'; end ode45set.refine = 8; ode45set.tol = dud.settings.tol; ode45set.solver = 'ode45'; ode45set.stepsize = dud.settings.stepsize; ode45set.hmax = 0; if strcmp(solver,'ode45') ode45ch = 'on'; else ode45ch = 'off'; end ode23tset.refine = 1; ode23tset.tol = dud.settings.tol; ode23tset.solver = 'ode23t'; ode23tset.stepsize = dud.settings.stepsize; ode23tset.hmax = 0; if strcmp(solver,'ode23t') ode23tch = 'on'; else ode23tch = 'off'; end ode23tbset.refine = 1; ode23tbset.tol = dud.settings.tol; ode23tbset.solver = 'ode23tb'; ode23tbset.stepsize = dud.settings.stepsize; ode23tbset.hmax = 0; if strcmp(solver,'ode23tb') ode23tbch = 'on'; else ode23tbch = 'off'; end dud.solver(1) = uimenu(mesolve,'label','Dormand Prince',... 'checked',dpch,... 'call',solverstr,'UserData',dpset); dud.solver(2) = uimenu(mesolve,'label','Runge-Kutta 4',... 'checked',rk4ch,... 'call',solverstr,'UserData',rk4set); dud.solver(3) = uimenu(mesolve,'label','ode45',... 'checked',ode45ch,... 'separator','on',... 'call',solverstr,'UserData',ode45set); dud.solver(4) = uimenu(mesolve,'label','ode23',... 'checked',ode23ch,... 'call',solverstr,'UserData',ode23set); dud.solver(5) = uimenu(mesolve,'label','ode113',... 'checked',ode113ch,... 'call',solverstr,'UserData',ode113set); dud.solver(6) = uimenu(mesolve,'label','ode15s',... 'checked',ode15sch,... 'call',solverstr,'UserData',ode15sset); dud.solver(7) = uimenu(mesolve,'label','ode23s',... 'checked',ode23sch,... 'call',solverstr,'UserData',ode23sset); dud.solver(8) = uimenu(mesolve,'label','ode23t',... 'checked',ode23tch,... 'call',solverstr,'UserData',ode23tset); dud.solver(9) = uimenu(mesolve,'label','ode23tb',... 'checked',ode23tbch,... 'call',solverstr,'UserData',ode23tbset); plotch = [ 'dud = get(gcf,''UserData'');',... 'aud = get(dud.axes,''UserData'');',... 'if aud.plot,',... ' aud.plot = 0;',... ' set(gcbo,''label'',''Plot while computing'');',... 'else,',... ' aud.plot = 1;',... ' set(gcbo,''label'',''Do not plot while computing'');',... 'end,',... 'set(dud.axes,''UserData'',aud);']; medir = uimenu(menu(1),'label','Solution direction.'); directionstr = ['ud = get(gcf,''UserData'');',... 'me = gcbo;',... 'ud.dir = get(me,''UserData'');',... 'set(ud.direction,''checked'',''off'');',... 'set(me,''checked'',''on'');',... 'set(gcf,''UserData'',ud);']; dud.direction(1) = uimenu(medir,'label','Both',... 'checked','on',... 'UserData',0,... 'call',directionstr); dud.dir = 0; dud.direction(2) = uimenu(medir,'label','Forward',... 'UserData',1,... 'call',directionstr); dud.direction(3) = uimenu(medir,'label','Back',... 'UserData',-1,... 'call',directionstr); markstr = ['ud = get(gcf,''UserData'');',... 'me = gcbo;',... 'chkd = get(me,''checked'');',... 'if strcmp(chkd,''on''),',... ' set(me,''checked'',''off'');',... ' ud.markflag = 0;',... 'else,',... ' set(me,''checked'',''on'');',... ' ud.markflag = 1;',... 'end,',... 'set(gcf,''UserData'',ud);']; dud.markflag = 0; dud.mark = uimenu(menu(1),'label','Mark initial points.',... 'checked','off',... 'call',markstr); meexportdata = uimenu(menu(1),'label','Export solution data.',... 'call','pplane81(''export'')',... 'separator','off','tag','dexp'); meplot = uimenu(menu(1),'label','Do not plot while computing',... 'call',plotch,'separator','on'); menu(6) = uimenu(menu(1),'label','Make the Display Window inactive.',... 'call','pplane81(''hotcold'')','separator','on'); dud.menu = menu; % View menu set(findobj(gcf,'label','&View'),'pos',6); set(findobj(gcf,'label','&Figure Toolbar'),... 'call','pplane81(''showbar'')'); set(0,'showhiddenhandles',hhsetup); set(gcf,'WindowButtonDownFcn','pplane81(''down'')'); set(ppdisp,'WindowButtonMotionFcn','pplane81(''cdisp'')'); hh1 = [dud.axes,dud.title.axes,nframe,dud.notice,dbutt([1 2 3])]; set(hh1,'units','norm'); hh2 = [nframe,dud.notice,dbutt(2:3),dud.axes,dud.menu(1),meplot,megraph]; set(hh2,'visible','on'); set(ppdisp,'vis','on'); dud.printstr = 'print -noui'; end % if (~isempty(ppdisp)) & else ppdispa = dud.axes; axes(ppdispa); cla xlabel(txname); ylabel(tyname); % Initialize the window matrix. set(findobj('tag','zbmenu'),'enable','off'); tstr1 = [txname,' '' = ', txderstr]; tstr2 = [tyname,' '' = ', tyderstr]; tstr = str2mat(tstr1,tstr2); dud.tstr = tstr; k = find(pflag); if ~isempty(k) lk = length(k); switch lk case 1 pstr1 = [pstring(k);{' '}]; case 2 pstr1 = pstring(k); case 3 pstr1 = pstring(k([2,3])); pstr2 = [pstring(k(1));{' '}]; case 4 pstr1 = pstring(k([2,4])); pstr2 = pstring(k([1,3])); case 5 pstr1 = pstring(k([3,5])); pstr2 = pstring(k([2,4])); pstr3 = [pstring(k(1));{' '}]; case 6 pstr1 = pstring(k([3,6])); pstr2 = pstring(k([2,5])); pstr3 = pstring(k([1,4])); end end set(dud.title.eq,'string',tstr); set(dud.title.p1,'string',pstr1); ext = get(dud.title.p1,'extent'); pos = get(dud.title.p1,'pos'); p1 = min(.9, .93 - ext(3)); pos(1) = p1; set(dud.title.p1,'pos',pos); set(dud.title.p2,'string',pstr2); ext = get(dud.title.p2,'extent'); pos = get(dud.title.p2,'pos'); p2 = min(.8, p1 - ext(3)-0.02); pos(1) = p2; set(dud.title.p2,'pos',pos); set(dud.title.p3,'string',pstr3); ext = get(dud.title.p3,'extent'); pos = get(dud.title.p3,'pos'); pos(1) = min(.7, p2 - ext(3)-0.02); set(dud.title.p3,'pos',pos); % Initialize important information as user data. dud.function = dfcn; dud.solhand = []; % Handles to solution curves. dud.ephand = []; % Handles to equilibrium points. dud.arr = []; % Handles for the direction and vector fields. dud.eqpts = []; % Equilibrium point data. dud.ics = []; % Marked initial points. dud.wmat = []; dud.color = sud.color; set(ppdisp,'UserData',dud); if strcmp(dud.syst.fieldtype,'nullclines') menunull = findobj('tag','null'); set(menunull,'label','Hide nullclines.'); end ud.y = zeros(2,1); ud.i = 0; ud.line = 0; wind = dud.syst.wind(:); dwind = [wind(1); wind(3); -wind(2); -wind(4)]; ud.DY = [wind(2)-wind(1); wind(4)-wind(3)]; ud.cwind = dwind - dud.settings.magn*[ud.DY;ud.DY]; ud.R = zeros(2,2); ud.rr = zeros(2,2); ud.perpeps = 0; ud.paraeps = 0; ud.sinkeps = 0; ud.minNsteps = 0; ud.turn = zeros(2,10); ud.tk = 0; ud.stop = 0; ud.gstop = 1; ud.plot = 1; set(dud.axes,'UserData',ud); tc = findobj('tag','testcase'); if ~isempty(tc) if strcmp(dud.syst.name,'default system') set(tc,'vis','on'); else set(tc,'vis','off'); end end ppkbd = findobj('name','pplane81 Keyboard input','vis','on'); if ~isempty(ppkbd),pplane81('kbd'),end pplevel = findobj('name','pplane81 Level sets'); if ~isempty(pplevel),delete(pplevel),end set(dud.title.axes,'handle','off'); pplane81('dirfield',ppdisp); end % if sud.flag == 1 & else elseif strcmp(action,'linear') % Initialize linearization window. % Linear takes the information from the EQPT window and initializes % the Linear Display window if it already exists. If the Linear % Display window does not exist, it builds one. ppeqpt = gcf; % Get the information from pplane81 Equilibrium ... sud = get(ppeqpt,'UserData'); WINvect = [-1 1 -1 1]; jac = sud.jac; type = sud.type; vectors = sud.vectors; ppdisp = findobj('name','pplane81 Display'); pdud = get(ppdisp,'UserData'); settings = pdud.settings; system = sud.system; % Find pplane81 Linearization if it exists. % If pplane81 Linearization exists, update it. If it does not build it. pplin= findobj('name','pplane81 Linearization'); if (~isempty(pplin)) figure(pplin); dud = get(pplin,'UserData'); dud.syst = system; dud.settings = settings; dfcn = dud.function; else pplin = figure('name','pplane81 Linearization',... 'NumberTitle','off',... 'interrupt','on',... 'visible','off',... 'tag','pplane81'); pplane81('figdefault',pplin); dud = get(pplin,'UserData'); dud.dir = 0; dud.syst = system; dud.settings = settings; dfcn = ''; fs = dud.fontsize; ssize = dud.ssize; r = ssize/10; ppaxw = 300; % Default axes width ppaxh = 300; % Default axes height ppaxl = 45; % Default axes left buttw = 40; % Default button width titleh = 45; % Default title height ppaxb = 35; ppdh = ppaxb + ppaxh + titleh; uni = get(0,'units'); set(0,'units','pixels'); ss = get(0,'screensize'); set(0,'units',uni); sw = ss(3);sh = ss(4); if r*ppdh > sh -80; r = (sh-80)/ppdh; fs = fs*r; lw = 0.5*r; set(gcf,'defaultaxesfontsize',fs,'defaultaxeslinewidth',lw); set(gcf,'defaulttextfontsize',fs); set(gcf,'defaultlinelinewidth',lw); set(gcf,'defaultuicontrolfontsize',fs*0.9); end axpos = r*[ppaxl,ppaxb,ppaxw,ppaxh]; % Set up the plot axes. dud.axes = axes('units','pix',... 'position',axpos,... 'next','add',... 'box','on',... 'interrupt','on',... 'xgrid','on',... 'ygrid','on',... 'drawmode','fast',... 'visible','off',... 'tag','linear axes'); % Set up the title. dud.title.axes = axes('box','off','xlim',[0 1],'ylim',[0 1],... 'units','pix','vis','off','xtick',[-1],'ytick',[-1],... 'xticklabel','','yticklabel',''); dud.title.eq = text(0.07,0.5,' ','vert','middle'); dud.title.p1 = text(0.75,0.5,' ','vert','middle'); dud.title.p2 = text(0.65,0.5,' ','vert','middle'); tstr = {'x_2';'y^2'}; set(dud.title.eq,'string',tstr,'units','pix'); ext = get(dud.title.eq ,'extent'); titleh = ext(4)+15*r; set(dud.title.eq,'units','data'); taxpos = r*[ppaxl,ppaxb+ppaxh,ppaxw,titleh]; set(dud.title.axes,'pos',taxpos,'color',get(gcf,'color')); % Finish the positions. % buttw = r*buttw; % butth = fs+10*r; % buttl = ppaxl + ppaxw + 5; ppdw = r*(ppaxl + ppaxw + 35); ppdh = r*(ppaxb+ppaxh)+titleh; ppdleft = sw-ppdw-60; ppdbot = 60; ppdpos = [ppdleft,ppdbot,ppdw,ppdh]; set(pplin,'resize','on'); set(pplin,'pos',ppdpos); Arrflag = system.fieldtype; axpos = r*[ppaxl,ppaxb,ppaxw,ppaxh]; stopstr = 'aud = get(gca,''UserData'');aud.stop = 4;set(gca,''UserData'',aud);'; stoppos = [axpos(1)+axpos(3)-r*10, axpos(2)+axpos(4)+r*10, r*40,fs+10*r]; dud.butt = uicontrol('style','push',... 'pos',stoppos,... 'string','Stop','call',stopstr,'vis','off','tag','stop'); hhsetup = get(0,'showhiddenhandles'); set(0,'showhiddenhandles','on'); % Configure the Toolbar. set(pplin,'ToolBar','none'); tmenu = findobj(gcf,'label','&Tools'); delete(tmenu); % File menu fmenu = findobj(gcf,'label','&File'); delete(findobj(fmenu,'label','&New Figure')); delete(findobj(fmenu,'label','&Open...')); delete(findobj(fmenu,'label','Pre&ferences...')); set(findobj(fmenu,'label','&Close'),'pos',1); set(findobj(fmenu,'label','&Save'),'pos',2,... 'separator','off'); set(findobj(fmenu,'label','Save &As...'),'pos',3); set(findobj(fmenu,'label','&Export...'),... 'pos',4); delete(findobj(fmenu,'label','Pre&ferences...')); set(findobj(fmenu,'label','Pa&ge Setup...'),'pos',5); set(findobj(fmenu,'label','Print Set&up...'),'pos',6); set(findobj(fmenu,'label','Print Pre&view...'),'pos',7); set(findobj(fmenu,'label','&Print...'),'pos',8); mequit = uimenu(fmenu,'label','Quit pplane81',... 'call','pplane81(''quit'')','separator','on','pos',9); % View menu set(findobj(gcf,'label','&Figure Toolbar'),... 'call','pplane81(''showbar'')'); % Solutions menu solmenu = uimenu('label','Solutions','pos',3); menukey = uimenu(solmenu,'label','Keyboard input.','call',... 'pplane81(''kbd'')','vis','on'); mesev = uimenu(solmenu,'label','Plot several solutions.',... 'call','pplane81(''several'')'); fundcall = ['dud = get(gcf,''UserData'');',... 'col = dud.color.sep;',... 'vect = dud.vectors;',... 'h = zeros(1,2);',... 'h(1) = plot(2*[vect(1,1),-vect(1,1)],2*[vect(2,1),-vect(2,1)]);',... 'h(2) = plot(2*[vect(1,2),-vect(1,2)],2*[vect(2,2),-vect(2,2)]);',... 'set(h,''color'',col);',... 'dud.solhand = [dud.solhand;h(:)];'... 'set(gcf,''UserData'',dud);']; dud.menu = uimenu(solmenu,'call',fundcall); markstr = ['ud = get(gcf,''UserData'');',... 'me = gcbo;',... 'chkd = get(me,''checked'');',... 'if strcmp(chkd,''on''),',... ' set(me,''checked'',''off'');',... ' ud.markflag = 0;',... 'else,',... ' set(me,''checked'',''on'');',... ' ud.markflag = 1;',... 'end,',... 'set(gcf,''UserData'',ud);']; dud.markflag = pdud.markflag; chkd = 'off'; if dud.markflag chkd = 'on'; end dud.mark = uimenu(solmenu,'label','Mark initial points.',... 'checked',chkd,... 'call',markstr); % Edit menu emenu = findobj(gcf,'label','&Edit'); medallsol = uimenu(emenu,'label','Erase all solutions.',... 'call','pplane81(''dallsol'')',... 'pos',1); medallics = uimenu(emenu,'label','Erase all marked initial points.',... 'call','pplane81(''dallics'')',... 'separator','off',... 'pos',2); medall = uimenu(emenu,'label','Erase all graphics objects.',... 'call','pplane81(''dall'')',... 'separator','off',... 'pos',3); medel = uimenu(emenu,'label','Delete a graphics object.',... 'call','pplane81(''delete'')',... 'visible','on',... 'pos',4); menutext = uimenu(emenu,'label','Enter text on the Display Window.',... 'call','pplane81(''text'')',... 'pos',5); set(findobj(emenu,'label','&Undo'),'separator','on',... 'pos',6); set(findobj(emenu,'label','Cu&t'),'pos',7); set(findobj(emenu,'label','&Copy'),'pos',8); set(findobj(emenu,'label','&Paste'),'pos',9); set(findobj(emenu,'label','Clea&r'),'pos',10); set(findobj(emenu,'label','&Select All'),'pos',11); set(findobj(emenu,'label','Copy &Figure'),'pos',12); set(findobj(emenu,'label','Copy &Options'),'pos',13); set(findobj(emenu,'label','F&igure Properties'),'pos',14); set(findobj(emenu,'label','&Axes Properties'),'pos',15); set(findobj(emenu,'label','C&urrent Object Properties'),'pos',16); % Options menu optmenu = uimenu('label','Options','visible','off'); mehc = uimenu(optmenu,'label','Make the Display Window inactive.',... 'call','pplane81(''hotcold'')','separator','on'); set(0,'showhiddenhandles',hhsetup); set(pplin,'WindowButtonDownFcn',['pplane81(''down'',',num2str(pplin),')']); hh1 = [dud.axes,dud.title.axes]; set(hh1,'units','norm'); set(dud.axes,'visible','on'); set(pplin,'vis','on'); dud.printstr = 'print'; end % if (~isempty(pplin)) & else dud.ics = []; switch type case 1 vstr = 'on'; lstr = 'Plot stable and unstable orbits.'; case {2,3} vstr = 'on'; lstr = 'Plot fundamental modes.'; otherwise vstr = 'off'; lstr=''; end set(dud.menu,'label',lstr,'vis',vstr); pplina = dud.axes; dud.vectors = sud.vectors; axes(pplina); cla xlabel('u'); ylabel('v'); % The title strings. tstr = {'u'' = A u + B v';'v'' = C u + D v'}; dud.tstr = tstr; pstr2 = {['A = ',num2str(jac(1,1))];['C = ',num2str(jac(2,1))]}; pstr1 = {['B = ',num2str(jac(1,2))];['D = ',num2str(jac(2,2))]}; set(dud.title.eq,'string',tstr); set(dud.title.p1,'string',pstr1); ext = get(dud.title.p1,'extent'); pos = get(dud.title.p1,'pos'); p1 = min(.9, .93 - ext(3)); pos(1) = p1; set(dud.title.p1,'pos',pos); set(dud.title.p2,'string',pstr2); ext = get(dud.title.p2,'extent'); pos = get(dud.title.p2,'pos'); pos(1) = min(.8, p1 - ext(3)-0.02); set(dud.title.p2,'pos',pos); if (~strcmp(dfcn,'') & exist(dfcn)==2) delete([tempdir,dfcn,'.m']);end tee = clock; tee = ceil(tee(6)*100); dfcn=['pptp',num2str(tee)]; fcnstr = ['function ypr = ',dfcn,'(t,y)\n\n']; commstr = '%%%% Created by pplane81\n\n'; astr = ['A = [',num2str(jac(1,1)),',',num2str(jac(1,2)),';',... num2str(jac(2,1)),',',num2str(jac(2,2)),'];\n']; dstr = 'ypr = A * y;'; ppf = fopen([tempdir,dfcn,'.m'],'w'); fprintf(ppf,fcnstr); fprintf(ppf,commstr); fprintf(ppf,astr); fprintf(ppf,dstr); fclose(ppf); % Initialize important information as user data. dud.function = dfcn; dud.solhand = []; % Handles to solution curves. dud.ephand = []; % Handles to equilibrium points. dud.contours = []; dud.arr = []; % Handles for the direction and vector % fields. dud.eqpts = []; % Equilibrium point data. dud.notice = 0; dud.syst.wind = [-1 1 -1 1]'; dud.color = sud.color; ud.y = zeros(2,1); ud.i = 0; ud.line = 0; wind = [-1 1 -1 1]'; dwind = [wind(1); wind(3); -wind(2); -wind(4)]; ud.DY = [wind(2)-wind(1); wind(4)-wind(3)]; ud.cwind = dwind - dud.settings.magn*[ud.DY;ud.DY]; ud.R = zeros(2,2); ud.rr = zeros(2,2); ud.perpeps = 0; ud.paraeps = 0; ud.sinkeps = 0; ud.turn = zeros(2,10); ud.tk = 0; ud.stop = 0; ud.gstop = 1; ud.plot = 1; dud.ephand = plot(0,0,... 'color',dud.color.eqpt,... 'markersize',20,... 'marker','.'); set(dud.axes,'UserData',ud); set(pplin,'UserData',dud); ppkbd = findobj('name','pplane81 Keyboard input','vis','on'); if ~isempty(ppkbd),pplane81('kbd'),end pplane81('dirfield',pplin); elseif strcmp(action,'dirfield') % 'dirfield' computes and plots the field elements. This is the entry % point both from 'display' and from later commands that require the % recomputation of the field elements. % Find pplane81 Display and get the user data. disph = input1; % This could be ppdisp or pplin. dud = get(disph,'UserData'); color = dud.color; dfcn = dud.function; ppdispa = dud.axes; WINvect = dud.syst.wind; settings = dud.settings; notice = dud.notice; if notice ~= 0 nstr = get(notice,'string'); nstr(1:4)=nstr(2:5); nstr{5,1} = 'Computing the field elements.'; set(notice,'string',nstr); % Augment the window matrix wmat = dud.wmat; wrows = size(wmat,1); wflag = 0; for k = 1:wrows if wmat(k,:)==WINvect wflag = 1; end end if wflag == 0 wmat = [wmat;WINvect]; dud.wmat = wmat; end if wrows hhsetup = get(0,'showhiddenhandles'); set(0,'showhiddenhandles','on'); hhh = findobj('tag','zbmenu'); set(hhh,'enable','on'); set(0,'showhiddenhandles',hhsetup); end end Xmin = WINvect(1); Xmax = WINvect(2); Ymin = WINvect(3); Ymax = WINvect(4); N = dud.syst.npts; if strcmp(get(disph,'name'),'pplane81 Linearization') N = floor(3*N/4); end deltax=(Xmax - Xmin)/(N-1); deltay=(Ymax - Ymin)/(N-1); % Set up the display window. Dxint=[Xmin-deltax,Xmax+deltax]; Dyint=[Ymin-deltay,Ymax+deltay]; % Set up the original mesh. XXXg=Xmin + deltax*[0:N-1]; YYYg=Ymin + deltay*[0:N-1]; [Xx,Yy]=meshgrid(XXXg,YYYg); % Calculate the line and vector fields. Xx=Xx(:);Yy=Yy(:); Ww = zeros(size(Xx)); Ww = feval(dfcn,0,[Xx';Yy']); Vv = Ww(1,:) + Ww(2,:)*sqrt(-1); Vv = Vv.'; Arrflag = dud.syst.fieldtype; mgrid = Xx+Yy.*sqrt(-1); % mgrid = mgrid(:); zz=Vv.'; sc = min(deltax,deltay); arrow=[-1,1].'; zzz=sign(zz); scale = sqrt((real(zzz)/deltax).^2+(imag(zzz)/deltay).^2); ww = (zzz == 0); scale = scale + ww; aa1 = 0.3*arrow*(zzz./scale)+ones(size(arrow))*(mgrid.'); [r,c] = size(aa1); aa1=[aa1;NaN*ones(1,c)]; aa1=aa1(:); arrow = [0,1,.7,1,.7].' + [0,0,.25,0,-.25].' * sqrt(-1); zz=sign(zz).*((abs(zz)).^(1/3)); scale = 0.9*sc./max(max(abs(zz))); aa2 = scale*arrow*zz +ones(size(arrow))*(mgrid.'); [r,c] = size(aa2); aa2=[aa2;NaN*ones(1,c)]; aa2=aa2(:); axes(ppdispa); arr = dud.arr; % Delete the old field data. if isstruct(arr) hand = [arr.hx;arr.hy;arr.lines;arr.arrows;arr.barrows]; delete(hand); end NN = N; N = 50; k = ceil(N/NN); deltax=(Xmax - Xmin)/(N-1); deltay=(Ymax - Ymin)/(N-1); % Set up the original mesh. XXXg=Xmin + deltax*[-k:N+k]; YYYg=Ymin + deltay*[-k:N+k]; [Xx,Yy]=meshgrid(XXXg,YYYg); % Calculate the line and vector fields. Xxx=Xx(:);Yyy=Yy(:); Ww = zeros(size(Xxx)); Ww = feval(dfcn,0,[Xxx';Yyy']); DX = Ww(1,:)'; DY = Ww(2,:)'; minx = min(DX); maxx = max(DX); miny = min(DY); maxy = max(DY); DR = Xxx.*DX + Yyy.*DY; DTheta = Xxx.*DY - Yyy.*DX; minr = min(DR); maxr = max(DR); minth = min(DTheta); maxth = max(DTheta); DX = reshape(DX,N+2*k+1,N+2*k+1); DY = reshape(DY,N+2*k+1,N+2*k+1); DR = reshape(DR,N+2*k+1,N+2*k+1); DTheta = reshape(DTheta,N+2*k+1,N+2*k+1); if minx < 0 & 0 < maxx [Cx,hx] = contour(Xx,Yy,DX,[0,0],'--'); set(hx,'visible','off','color',color.xcline,'linestyle','--'); else hx = zeros(0,1); end if miny < 0 & 0 < maxy [Cy,hy] = contour(Xx,Yy,DY,[0,0],'--'); set(hy,'visible','off','color',color.ycline,'linestyle','--'); else hy = zeros(0,1); end if minr < 0 & 0 < maxr [Cr,hr] = contour(Xx,Yy,DR,[0,0],'--'); set(hr,'visible','off','color',color.xcline,'linestyle','--'); else hr = zeros(0,1); end if minth < 0 & 0 < maxth [Cth,hth] = contour(Xx,Yy,DTheta,[0,0],'--'); set(hth,'visible','off','color',color.ycline,'linestyle','--'); else hth = zeros(0,1); end arrh1 = plot(real(aa1),imag(aa1),'color',color.arrows,'visible','off'); arrh2 = plot(real(aa2),imag(aa2),'color',color.arrows,'visible','off'); % We plot both the line field and the vector field. Then we % control which is seen by manipulating the visibility. Xmin = Xmin - k*deltax; Xmax = Xmax + k*deltax; Ymin = Ymin - k*deltay; Ymax = Ymax + k*deltay; Zz = (Xx-Xmin).*(Xmax-Xx).*(Yy-Ymin).*(Ymax-Yy).*abs(DX).*abs(DY); Zzc = Zz(2:N+k+2,2:N+k+2); Xxc = Xx(2:N+k+2,2:N+k+2); Yyc = Yy(2:N+k+2,2:N+k+2); DXc = DX(2:N+k+2,2:N+k+2); DYc = DY(2:N+k+2,2:N+k+2); Zzl = Zz(1:N+k+1,2:N+k+2); Zzr = Zz(3:N+k+3,2:N+k+2); Zzd = Zz(2:N+k+2,1:N+k+1); Zzu = Zz(2:N+k+2,3:N+k+3); Zzdl = Zz(1:N+k+1,1:N+k+1); Zzdr = Zz(3:N+k+3,1:N+k+1); Zzul = Zz(1:N+k+1,3:N+k+3); Zzur = Zz(3:N+k+3,3:N+k+3); kk = find((Zzc>Zzl)&(Zzc>Zzr)&(Zzc>Zzd)&(Zzc>Zzu)... &(Zzc>Zzdl)&(Zzc>Zzdr)&(Zzc>Zzul)&(Zzc>Zzur)); Xxx = Xxc(kk); Yyy = Yyc(kk); DXx = DXc(kk); DYy = DYc(kk); ll = length(hx); Xxx = Xxx(:); Yyy = Yyy(:); DXx = DXx(:); DYy = DYy(:); if ~exist('both') Xxx = zeros(0,1); Yyy = zeros(0,1); DXx = zeros(0,1); DYy = zeros(0,1); end for j = 1:ll xd = get(hx(j),'xdata');xd = xd(:); yd = get(hx(j),'ydata');yd = yd(:); zxd = feval(dfcn,0,[xd';yd']); yxd = zxd(2,:);yxd = yxd(:); ayxd = abs(yxd).*(xd-Xmin).*(Xmax-xd).*(yd-Ymin).*(Ymax-yd); NNN = length(ayxd); ayxdc = ayxd(2:NNN-1); yxdc = yxd(2:NNN-1); ayxdl = ayxd(1:NNN-2); ayxdr = ayxd(3:NNN); lll = find((ayxdc>ayxdl)&(ayxdc>ayxdr)); Xx = xd(lll+1); Yy = yd(lll+1); Xxx = [Xxx;Xx(:)]; Yyy = [Yyy;Yy(:)]; yxdcp = yxdc(lll); yxdcp = yxdcp(:); DXx = [DXx;zeros(size(yxdcp))]; DYy = [DYy;yxdcp]; end ll = length(hy); for j = 1:ll xd = get(hy(j),'xdata');xd = xd(:); yd = get(hy(j),'ydata');yd = yd(:); zxd = feval(dfcn,0,[xd';yd']); yxd = zxd(1,:);yxd = yxd(:); ayxd = abs(yxd).*(xd-Xmin).*(Xmax-xd).*(yd-Ymin).*(Ymax-yd); NNN = length(ayxd); ayxdc = ayxd(2:NNN-1); yxdc = yxd(2:NNN-1); ayxdl = ayxd(1:NNN-2); ayxdr = ayxd(3:NNN); lll = find((ayxdc>ayxdl)&(ayxdc>ayxdr)); Xx = xd(lll+1); Yy = yd(lll+1); Xxx = [Xxx;Xx(:)]; Yyy = [Yyy;Yy(:)]; yxdcp = yxdc(lll); yxdcp = yxdcp(:); DXx = [DXx;yxdcp]; DYy = [DYy;zeros(size(yxdcp))]; end mgrid = Xxx(:) + Yyy(:)*sqrt(-1); Vv = DXx(:) + DYy(:)*sqrt(-1); Vv = sign(Vv); %.*((abs(Vv)).^(1/3)); sc = (N/20)*min(deltax,deltay); aa3 = sc*arrow*(Vv.') +ones(size(arrow))*(mgrid.'); [r,c] = size(aa3); aa3 = [aa3;NaN*ones(1,c)]; aa3=aa3(:); arrh3 = plot(real(aa3),imag(aa3),'color',color.narrows,'visible','off'); switch Arrflag case 'nullclines' % set([hx;hy;arrh3],'vis','on'); set([hx;hy],'vis','on'); case 'lines' set(arrh1,'visible','on'); case 'arrows' set(arrh2,'visible','on'); end dud.arr.lines = arrh1; % Save the handles for later use. dud.arr.arrows = arrh2; % Save the handles for later use. dud.arr.hx = hx; dud.arr.hy = hy; % dud.arr.barrows = arrh3; dud.arr.barrows = []; dud.arr.hr = hr; dud.arr.hth = hth; menull = findobj('tag','null'); mernull = findobj('tag','rnull'); if strcmp(get(menull,'label'),'Delete nullclines.') set([hx;hy],'vis','on'); elseif strcmp(get(mernull,'label'),'Delete polar nullclines.') set([hr;hth],'vis','on'); end if notice ~= 0 nstr = get(notice,'string'); nstr(1:4) = nstr(2:5); nstr{5,1} = 'Ready.'; set(notice,'string',nstr); end set(disph,'UserData',dud); axis([Dxint,Dyint]); elseif strcmp(action,'hotcold') % 'hotcold' is the callback for the menu selection that makes the % Display Window active or inactive. ppdisp = gcf; dud = get(ppdisp,'UserData'); nstr = get(dud.notice,'string'); nstr(1:4) = nstr(2:5); mehc = dud.menu(6); if (findstr(get(mehc,'label'),'inactive')) set(ppdisp,'WindowButtonDownFcn',' '); set(mehc,'label','Make the Display Window active.'); nstr{5,1} = 'The Display Window is not active.'; set(dud.notice,'string',nstr); else set(ppdisp,'WindowButtonDownFcn','pplane81(''down'')'); set(mehc,'label','Make the Display Window inactive.'); nstr{5,1} = 'The Display Window is active.'; set(dud.notice,'string',nstr); end elseif strcmp(action,'down') % 'down' is the Window Button Down call. It starts the computation of % solutions from a click of the mouse. disph = gcf; seltype = get(disph,'selectiontype'); if strcmp(seltype,'alt') pplane81('zoom'); return elseif strcmp(seltype,'extend') pplane81('zoomsqd'); return end dud = get(disph,'UserData'); ax = dud.axes; ch = findobj('type','uicontrol','enable','on'); set(ch,'enable','inactive'); wbdf = get(disph,'WindowbuttonDownFcn'); set(disph,'WindowbuttonDownFcn',''); axes(ax); initpt = get(ax,'currentpoint'); initpt = initpt(1,[1,2]); if dud.markflag h = plot(initpt(1),initpt(2),'.k'); dud.ics = [dud.ics,h]; set(disph,'UserData',dud); end pplane81('solution',initpt,disph); set(disph,'WindowbuttonDownFcn',wbdf); % set([ch;mh],'enable','on'); set(ch,'enable','on'); notice = dud.notice; if notice ~= 0 nstr = get(notice,'string'); nstr(1:4) = nstr(2:5); nstr{5,1} = 'Ready.'; set(notice,'string',nstr) end elseif strcmp(action,'several') % 'several' allows the user to pick several initial points at once. % This is not needed in X-windows, but it is on the Macintosh. disph = gcf; ch = findobj('type','uicontrol','enable','on'); % mh = findobj('type','uimenu','enable','on'); set(ch,'enable','inactive'); % set(mh,'enable','off') wbdf = get(disph,'WindowbuttonDownFcn'); set(disph,'WindowbuttonDownFcn',''); dud = get(disph,'UserData'); notice = dud.notice; if notice ~= 0 nstr = get(notice,'string'); nstr(1:4) = nstr(2:5); nstr{5,1} = 'Pick initial points with the mouse. Enter "Return" when finished.'; set(notice,'string',nstr) end % [X,Y]=ppginput; [X,Y]=ginput; NN = length(X); for k = 1:NN initpt = [X(k),Y(k)]; pplane81('solution',initpt,disph); end if notice ~= 0 nstr = get(notice,'string'); nstr(1:4) = nstr(2:5); nstr{5,1} = 'Ready.'; set(notice,'string',nstr); end set(ch,'enable','on'); % set([ch;mh],'enable','on'); set(disph,'WindowbuttonDownFcn',wbdf); elseif strcmp(action,'test case') tic ppdisp = gcf; for k = 1:10 initpt = k*[-.2 .2]; pplane81('solution',initpt,ppdisp); initpt = k*[-.2 -.2]; pplane81('solution',initpt,ppdisp); end pplane81('solution',[1,-1],ppdisp); dud = get(gcf,'UserData'); nstr = get(dud.notice,'string'); nstr(1:4) = nstr(2:5); nstr{5,1} = 'Ready.'; set(dud.notice,'string',nstr) toc elseif strcmp(action,'solution') % 'solution' effects the computation and (erasemode == none) plotting of % solutions. It also stores the data as appropriate. disph = input2; dud = get(disph,'UserData'); tcol = dud.color.temp; pcol = dud.color.orb; notice = dud.notice; initpt = input1(:); dfcn = dud.function; ppdispa = dud.axes; settings = dud.settings; ptstr = [' (',num2str(initpt(1),2), ', ', num2str(initpt(2),2), ')']; refine = settings.refine; ssize = settings.stepsize; tol = settings.tol; ud = get(dud.axes,'UserData'); % rtol = tol; atol = tol*ud.DY*1e-4'; if length(initpt) == 2 AA = -1e6; BB = 1e6; switch dud.dir case 0 intplus = [0, BB]; intminus = [0, AA]; case -1 intplus = [0, 0]; intminus = [0, AA]; case 1 intplus = [0, BB]; intminus = [0, 0]; end else intplus = [initpt(3),initpt(5)]; intminus = [initpt(3),initpt(4)]; initpt = initpt([1:2]); end stopbutt = findobj(disph,'tag','stop'); set(stopbutt,'vis','on','enable','on'); solver = settings.solver; switch solver case 'Dormand Prince' solh = @ppdp45; opt = disph; case 'Runge-Kutta 4' solh = @pprk4; opt = disph; case 'ode45' solh = @ode45; opt = odeset('OutputFcn',@ppout,'Refine',refine,... 'RelTol',tol,'Abstol',atol); case 'ode23' solh = @ode23; opt = odeset('OutputFcn',@ppout,'Refine',refine,... 'RelTol',tol,'Abstol',atol); case 'ode113' solh = @ode113; opt = odeset('OutputFcn',@ppout,'Refine',refine,... 'RelTol',tol,'Abstol',atol); case 'ode15s' solh = @ode15s; opt = odeset('OutputFcn',@ppout,'Refine',refine,... 'RelTol',tol,'Abstol',atol); case 'ode23s' solh = @ode23s; opt = odeset('OutputFcn',@ppout,'Refine',refine,... 'RelTol',tol,'Abstol',atol); case 'ode23t' solh = @ode23t; opt = odeset('OutputFcn',@ppout,'Refine',refine,... 'RelTol',tol,'Abstol',atol); case 'ode23tb' solh = @ode23tb; opt = odeset('OutputFcn',@ppout,'Refine',refine,... 'RelTol',tol,'Abstol',atol); end exist(dfcn); dfh = str2func(dfcn); cflag = 0; if intplus(2)>intplus(1) cflag = cflag + 1; if notice ~= 0 nstr = get(notice,'string'); nstr(1:4) = nstr(2:5); nstr{5} = ['The forward orbit from',ptstr]; set(notice,'string',nstr); end drawnow [tp,xp] = feval(solh,dfh,intplus,initpt,opt); aud = get(ppdispa,'UserData'); hnew1 = aud.line; end if intminus(2) < intminus(1) cflag = cflag + 2; if notice ~= 0 nstr = get(notice,'string'); nstr(1:4) = nstr(2:5); nstr{5} = ['The backward orbit from',ptstr]; set(notice,'string',nstr); end drawnow [tm,xm] = feval(solh,dfh,intminus,initpt,opt); aud = get(ppdispa,'UserData'); hnew2 = aud.line; set(stopbutt,'vis','off'); end % if intminus(2) < intminus(1) % Store the trajectory. switch cflag case 1 % positive only set(hnew1,'xdata',xp(:,1),'ydata',xp(:,2),'zdata',tp,'color',pcol); set(hnew1,'erase','normal'); dud.solhand = [dud.solhand;hnew1]; case 2 % negative only x = flipud(xm); t = flipud(tm); set(hnew2,'xdata',x(:,1),'ydata',x(:,2),'zdata',t,'color',pcol); set(hnew2,'erase','normal'); dud.solhand = [dud.solhand;hnew2]; case 3 % both directions x = flipud(xm); t = flipud(tm); x=[x;xp]; t=[t;tp]; delete(hnew1); set(hnew2,'xdata',x(:,1),'ydata',x(:,2),'zdata',t,'color',pcol); set(hnew2,'erase','normal'); dud.solhand = [dud.solhand;hnew2]; end % switch cflag set(disph,'UserData',dud); elseif strcmp(action,'kcompute') % 'kcompute' is the call back for the Compute % button on the pplane81 Keyboard figure. compute = 1; kh = get(gcf,'UserData'); ppdisp = kh.fig; if (isempty(ppdisp)) pplane81('confused'); end dud = get(ppdisp,'UserData'); ppdispa = dud.axes; aud = get(ppdispa,'UserData'); ppset = findobj('name','pplane81 Setup'); sud = get(ppset,'UserData'); ch = findobj('type','uicontrol','enable','on'); set(ch,'enable','inactive'); set(ppdisp,'WindowbuttonDownFcn',''); xstr = get(kh.xval,'string'); ystr = get(kh.yval,'string'); pnameh = sud.h.pname; pvalh = sud.h.pval; pflag = zeros(1,4); perr = []; for kk = 1:6; pn = char(get(pnameh(kk),'string')); pv = char(get(pvalh(kk),'string')); if ~isempty(pn) if isempty(pv) perr = pvalh(kk); else xstr = pplane81('paraeval',pn,pv,xstr); ystr = pplane81('paraeval',pn,pv,ystr); end end end xvalue = str2num(xstr); yvalue = str2num(ystr); if get(kh.spec,'value') tzero = str2num(get(kh.tval,'string')); t0 = str2num(get(kh.t0,'string')); tf = str2num(get(kh.tf,'string')); initpt = [xvalue,yvalue,tzero,t0,tf]; str1 = 'Values must be entered for all of the entries.'; if (length(initpt) ~= 5) warndlg({str1},'Illegal input'); compute = 0; elseif tf <= t0 warndlg({'The final time of the computation interval';... 'must be smaller than the initial time.'},'Illegal input'); compute = 0; elseif (tzero < t0) | (tzero > tf) str2 = 'The initial time must be in the computation interval.'; warndlg(str2,'Illegal input'); compute = 0; end aud.gstop = 0; set(ppdispa,'UserData',aud); else initpt = [xvalue,yvalue]; if (length(initpt) ~= 2) warndlg({str1},'Illegal input'); compute = 0; end end % if get(kh.spec,'value') if compute if dud.markflag figure(ppdisp) h = plot(initpt(1),initpt(2),'.k'); dud.ics = [dud.ics,h]; set(ppdisp,'UserData',dud); end pplane81('solution',initpt,ppdisp); end if dud.notice ~= 0 nstr = get(dud.notice,'string'); nstr(1:4) = nstr(2:5); nstr{5} = 'Ready.'; set(dud.notice,'string',nstr); end % set([ch;mh],'enable','on'); set(ch,'enable','on'); set(ppdisp,'WindowbuttonDownFcn','pplane81(''down'')'); aud.gstop = 1; set(ppdispa,'UserData',aud); elseif strcmp(action,'kbd') % 'kbd' is the callback for the Keyboard Input menu selection. It % sets up the pplane81 Keyboard figure which allows accurate input of % initial values using the keyboard. ppdisp = gcf; % The figure to be plotted in. dud = get(ppdisp,'UserData'); Xname = dud.syst.xvar; Yname = dud.syst.yvar; xnstr = ['The initial value of ',Xname,' = ']; ynstr = ['The initial value of ',Yname,' = ']; tnstr = 'The initial value of t = '; ppkbd = findobj('name','pplane81 Keyboard input'); if ~isempty(ppkbd) delete(ppkbd); end ppkbd = figure('name','pplane81 Keyboard input',... 'vis','off',... 'NumberTitle','off','tag','pplane81'); pplane81('figdefault',ppkbd); set(ppkbd,'menubar','none'); kbd.fr1 = uicontrol('style','frame'); kbd.fr2 = uicontrol('style','frame'); kbd.fr3 = uicontrol('style','frame'); dname = get(ppdisp,'name'); kbd.which = uicontrol('style','text','horiz','center',... 'string',['Data for ',dname]); kbd.inst = uicontrol('style','text','horiz','left',... 'string','Enter the initial conditions:'); kbd.xname = uicontrol('style','text',... 'horiz','right','string',xnstr); kbd.xval = uicontrol('style','edit',... 'string','','call',''); kbd.yname = uicontrol('style','text',... 'horiz','right',... 'string',ynstr); kbd.yval = uicontrol('style','edit',... 'string',''); kbd.tname = uicontrol('style','text',... 'horiz','right',... 'string',tnstr); kbd.tval = uicontrol('style','edit'); kbd.t0 = uicontrol('style','edit'); kbd.tf = uicontrol('style','edit'); kbd.t = uicontrol('style','text','string','<= t<= '); kbd.spec = uicontrol('style','check','horiz','center',... 'string','Specify a computation interval.'); kbd.comp = uicontrol('style','push',... 'string','Compute','call','pplane81(''kcompute'')'); kbd.close = uicontrol('style','push',... 'string','Close','call','set(gcf,''vis'',''off'')'); kbd.fig = ppdisp; left = 5; varl = 70; % buttw = 60; frsep = 1; nudge = 3; xex = get(kbd.xname,'extent'); ht = xex(4)+nudge; yex = get(kbd.yname,'extent'); nwdth = max(xex(3),yex(3)) + nudge; varl = varl*ht/19; fr1bot = 2*left + ht; fr1ht = 4*nudge + 3*ht; frw = 2*nudge + nwdth + varl; fr1wind = [left,fr1bot,frw,fr1ht]; set(kbd.fr1,'pos',fr1wind); tnb = fr1bot + nudge; nl = left+nudge; vl = nl + nwdth; tvwind = [vl,tnb,varl,ht]; tnwind = [nl,tnb,nwdth,ht]; set(kbd.tname,'pos',tnwind); set(kbd.tval,'pos',tvwind); intext = get(kbd.t,'extent'); tw = intext(3); margin = (frw -tw - 2*varl)/2; t0l = left+margin; intbot = tnb +ht + nudge; t0wind = [t0l,intbot,varl,ht]; set(kbd.t0,'pos',t0wind); tl = t0l + varl; tfl = tl + tw; twind = [tl,intbot,tw,ht]; tfwind = [tfl,intbot,varl,ht]; set(kbd.t,'pos',twind); set(kbd.tf,'pos',tfwind); specb = intbot + ht + nudge; specw = frw - 2*nudge; specwind = [nl,specb,specw,ht]; set(kbd.spec,'pos',specwind); fr2bot = fr1bot + fr1ht +frsep; ynb = fr2bot + nudge; xnb = ynb + ht; xnwind = [nl,xnb,nwdth,ht]; ynwind = [nl,ynb,nwdth,ht]; xvwind = [vl,xnb,varl,ht]; yvwind = [vl,ynb,varl,ht]; instb = xnb + ht + nudge; instw = nwdth + varl; instwind = [nl,instb,instw,ht]; fr2ht = 4*nudge + 3*ht; fr3bot = fr2bot + fr2ht+frsep; fr3ht = 2*nudge + ht; frw = 2*nudge + nwdth + varl; whichbot = fr3bot + nudge; whichw = frw - 2*nudge; whichwind = [nl,whichbot,whichw,ht]; fr3wind = [left,fr3bot,frw,fr3ht]; fr2wind = [left,fr2bot,frw,fr2ht]; figw = 2*left + frw; fight = 3*left + ht + fr1ht + fr2ht + fr3ht; figwind = [30,300,figw,fight]; buttw = (frw-left)/2; closel = left; compl = 2*left+buttw; clwind = [closel,left,buttw,ht]; compwind = [compl,left,buttw,ht]; set(ppkbd,'pos',figwind); set(kbd.fr2,'pos',fr2wind); set(kbd.fr3,'pos',fr3wind); set(kbd.which,'pos',whichwind); set(kbd.inst,'pos',instwind); set(kbd.xname,'pos',xnwind); set(kbd.yname,'pos',ynwind); set(kbd.xval,'pos',xvwind); set(kbd.yval,'pos',yvwind); set(kbd.comp,'pos',compwind); set(kbd.close,'pos',clwind); speccall = [ 'ud = get(gcf,''UserData'');',... 'if get(gcbo,''value''),',... ' set([ud.t0,ud.t,ud.tf,ud.tname,ud.tval],''enable'',''on'');',... 'else,',... ' set([ud.t0,ud.t,ud.tf,ud.tname,ud.tval],''enable'',''off'');',... 'end']; set(kbd.spec,'call',speccall); set(ppkbd,'resize','on'); set(findobj(ppkbd,'type','uicontrol'),'units','normal'); set([kbd.tval,kbd.t0],'string','0'); set(kbd.spec,'value',0); set(ppkbd,'UserData',kbd,'vis','on'); set([kbd.t0,kbd.t,kbd.tf,kbd.tname,kbd.tval],'enable','off') set(findobj(ppkbd,'type','uicontrol'),'units','normal'); edith = findobj(ppkbd,'style','edit'); set(edith,'backgroundcolor','w'); figure(ppkbd) elseif strcmp(action,'eqpt') % Find and classify equilibrium points. ppdisp = findobj('name','pplane81 Display'); dud = get(ppdisp,'UserData'); col = dud.color.eqpt; dfcn = dud.function; dbutt = dud.butt; menu = dud.menu; ppdispa =dud.axes; Dx = get(ppdispa,'xlim'); Dy = get(ppdispa,'ylim'); DY = [Dx(2)-Dx(1);Dy(2)-Dy(1)]; epsilon = 1e-4; nstr = get(dud.notice,'string'); nstr(1:4) = nstr(2:5); nstr{5} = 'Choose an approximation with the mouse.'; set(dud.notice,'string',nstr); % z0 = ppginput(1); z0 = ginput(1); z = pplane81('newton',z0,dfcn); flag = z(:,4); A = z(:,2:3); B = real(A); k = find(abs(B)<1e-8); B(k) = zeros(size(k)); C = imag(A); k = find(abs(C)<1e-8); C(k) = zeros(size(k)); A = B+C*sqrt(-1); z = z(:,1); if (~flag | norm((z-z0')./DY) > 1/5) nstr(1:4) = nstr(2:5); nstr{5} = ['There is not an equilibrium point near (',... num2str(z0(1)),', ', num2str(z0(2)), '). Ready']; set(dud.notice,'string',nstr); return end zero =find(abs(z) < epsilon); z(zero) = zeros(size(zero)); D=det(A); T=trace(A); if (D<-epsilon) string1 = ['There is a saddle point at (',... num2str(z(1)),', ', num2str(z(2)), ').']; string2 = ''; type = 1; elseif(D>epsilon) if (T*T>4*D+epsilon) if (T<0) string1 = ['There is a nodal sink at (',... num2str(z(1)),', ', num2str(z(2)), ').']; string2 = ''; type = 2; else string1 = ['There is a nodal source at (',... num2str(z(1)),', ', num2str(z(2)), ').']; string2 = ''; type = 3; end elseif(T*T<4*D-epsilon) if(T<-epsilon) string1 = ['There is a spiral sink at (',... num2str(z(1)),', ', num2str(z(2)), ').']; string2 = ''; type = 4; elseif(T>epsilon) string1 = ['There is a spiral source at (',... num2str(z(1)),', ', num2str(z(2)), ').']; string2 = ''; type = 5; else string1 = ['There is a spiral equilibrium point at (',... num2str(z(1)),', ', num2str(z(2)), ').']; string2 = 'Its specific type has not been determined.'; type = 6; end else if (T>epsilon) string1 = ['There is a source at (',... num2str(z(1)),', ', num2str(z(2)), ').']; string2 = 'Its specific type has not been determined.'; type = 7; elseif (T<-epsilon) string1 = ['There is a sink at (',... num2str(z(1)),', ', num2str(z(2)), ').']; string2 = 'Its specific type has not been determined.'; type = 8; else string1 = ['There is an equilibrium point at (',... num2str(z(1)),', ', num2str(z(2)), ').']; string2 = 'Its specific type has not been determined.'; type = 9; end end else string1 = ['There is an equilibrium point at (',... num2str(z(1)),', ', num2str(z(2)), ')']; string2 = 'Its specific type has not been determined.'; type = 9; end [V,D] = eig(A); EqPtList = dud.eqpts; infostr{1,1} = string1; infostr{2,1} = string2; infostr{3,1} = ' '; infostr{4,1} = 'The Jacobian is:'; infostr{5,1} = [' ',num2str(A(1,1)),' ',num2str(A(1,2))]; infostr{6,1} = [' ',num2str(A(2,1)),' ',num2str(A(2,2))]; infostr{7,1} = ''; infostr{8,1} = 'The eigenvalues and eigenvectors are:'; infostr{9,1} = [' ',num2str(D(1,1)),' (',num2str(V(1,1)),', ',num2str(V(2,1)),')']; infostr{10,1} = [' ',num2str(D(2,2)),' (',num2str(V(1,2)),', ',num2str(V(2,2)),')']; k=1; l=size(EqPtList,1); while ((k <= l)) if (norm((EqPtList(k,1:2)-z')./DY') <= 1e-3) break; end k = k+1; end if (k > l) EqPtList = [EqPtList;z',type]; dud.eqpts = EqPtList; newh = plot(z(1),z(2),... 'color',col,... 'markersize',20,... 'marker','.','Erasemode','none'); dud.ephand = [dud.ephand;newh]; set(newh,'Erasemode','normal'); drawnow end ppeqpt = findobj('name','pplane81 Equilibrium point data'); if (isempty(ppeqpt)) ppeqpt = figure('vis','off','resize','on',... 'name','pplane81 Equilibrium point data',... 'NumberTitle','off','tag','pplane81'); pplane81('figdefault',ppeqpt); set(ppeqpt,'menubar','none'); ud.frame = uicontrol('style','frame'); ud.eptext = uicontrol('style','text','string',infostr,'hor','left'); ud.goaway = uicontrol('style','push','string','Go away',... 'call','close'); ud.display = uicontrol('style','push',... 'string','Display the linearization',... 'call','pplane81(''linear'')'); else figure(ppeqpt); ud = get(ppeqpt,'UserData'); set(ppeqpt,'vis','off'); set(findobj(ppeqpt,'type','uicontrol'),'units','pixels'); set(ud.eptext,'string',infostr); end ud.jac = A; ud.vectors = V; ud.type = type; ud.system = dud.syst; ud.settings = dud.settings; ud.color = dud.color; set(ppeqpt,'UserData',ud); left = 5; nudge = 3; ext = get(ud.eptext,'extent'); n = size(infostr,1); ht = ext(4)/n+2*nudge; frbot = 3*left + 2*ht; txtbot = frbot + nudge; txtl = left + nudge; twind = [txtl,txtbot,ext(3)+nudge,ext(4)]; frw = ext(3) + 3*nudge; frh = ext(4) + 2*nudge; frwind = [left,frbot,frw,frh]; figw = frw + 2*left; figh = frh + 4*left + 2*ht; uni = get(0,'units'); set(0,'units','pixels'); ss = get(0,'screensize'); set(0,'units',uni); sh = ss(4); figbot = sh - figh - 350; figwind = [30,figbot,figw,figh]; buttw = frw - 4; buttl = left + 2; closewind = [buttl,left,buttw,ht]; dispb = 2*left + ht; dispwind = [buttl,dispb,buttw,ht]; set(ppeqpt,'pos',figwind); set(ud.frame,'pos',frwind); set(ud.eptext,'pos',twind); set(ud.goaway,'pos',closewind); set(ud.display,'pos',dispwind); set(findobj(ppeqpt,'type','uicontrol'),'units','normal'); set(ppeqpt,'vis','on'); set(get(ppeqpt,'child'),'vis','on'); nstr(1:4) = nstr(2:5); nstr{5} = 'Ready.'; set(dud.notice,'string',nstr); set(ppdisp,'UserData',dud); elseif strcmp(action,'stunst') ppdisp = findobj('name','pplane81 Display'); dud = get(ppdisp,'UserData'); ecol = dud.color.eqpt; scol = dud.color.sep; ppdispa = dud.axes; aud = get(ppdispa,'UserData'); DY = aud.DY; settings = dud.settings; dfcn = dud.function; EqPtList = dud.eqpts; Stop = norm(DY)*1e-4; % Plot the stable and unstable orbits at a saddle point. nstr = get(dud.notice,'string'); nstr(1:4) = nstr(2:5); nstr{5} = 'Choose a saddle point with the mouse.'; set(dud.notice,'string',nstr); % z0 = ppginput(1); z0 = ginput(1); z = zeros(2,1); z = pplane81('newton',z0,dfcn); flag = z(:,4); A = z(:,2:3); k = find(abs(A)<1e-8); A(k) = zeros(size(k)); z = z(:,1); if (~flag | norm((z-z0')./DY)> 1/5) nstr(1:4) = nstr(2:5); nstr{5} = ['There is not an equilibrium point near (',... num2str(z0(1),2),', ', num2str(z0(2),2), ').']; set(dud.notice,'string',nstr); return end zero = find(abs(z) < 1e-4); z(zero) = zeros(size(zero)); D=det(A); if (D>=0) nstr(1:4) = nstr(2:5); nstr{5} = ['The equilibrium point at (',... num2str(z(1),2),', ', num2str(z(2),2),... ') is not a saddle point.']; set(dud.notice,'string',nstr); return end if isempty(EqPtList) EqPtList = [z',1]; dud.eqpts = EqPtList; else k=1; l=size(EqPtList,1); while ((k <= l)) if (norm(EqPtList(k,1:2)-z') <= Stop) break; end k = k+1; end if (k > l) EqPtList = [EqPtList;z',1]; dud.eqpts = EqPtList; end end nstr(1:4) = nstr(2:5); nstr{5} = ['Plotting from the saddle point at (',... num2str(z(1),2),', ', num2str(z(2),2),').']; set(dud.notice,'string',nstr); [V,L] = eig(A); [L,I] = sort(diag(L)); magn = settings.magn; refine = settings.refine; ssize = settings.stepsize; solver = settings.solver; tol = settings.tol; offset = norm(DY)/1000; lm=abs(L); lm=lm/max(lm); lm = (max([lm,0.2*ones(size(lm))]'))'; offs = offset./lm; ud = get(dud.axes,'UserData'); atol = tol*ud.DY*1e-4'; options = odeset('OutputFcn',@ppout,'Refine',refine,... 'RelTol',tol,'Abstol',atol); stopbutt = findobj('tag','stop'); set(stopbutt,'vis','on','enable','on'); VV=V(:,I(1)); newhand = zeros(4,1); w = z + offs(1)*VV; nstr = get(dud.notice,'string'); nstr(1:4) = nstr(2:5); nstr{5} = 'The first stable trajectory'; set(dud.notice,'string',nstr); int = [0,-1e6]; solver = settings.solver; switch solver case 'Dormand Prince' solh = @ppdp45; opt = ppdisp; case 'Runge-Kutta 4' solh = @pprk4; opt = ppdisp; case 'ode45' solh = @ode45; opt = odeset('OutputFcn',@ppout,'Refine',refine,... 'RelTol',tol,'Abstol',atol); case 'ode23' solh = @ode23; opt = odeset('OutputFcn',@ppout,'Refine',refine,... 'RelTol',tol,'Abstol',atol); case 'ode113' solh = @ode113; opt = odeset('OutputFcn',@ppout,'Refine',refine,... 'RelTol',tol,'Abstol',atol); case 'ode15s' solh = @ode15s; opt = odeset('OutputFcn',@ppout,'Refine',refine,... 'RelTol',tol,'Abstol',atol); case 'ode23s' solh = @ode23s; opt = odeset('OutputFcn',@ppout,'Refine',refine,... 'RelTol',tol,'Abstol',atol); case 'ode23t' solh = @ode23t; opt = odeset('OutputFcn',@ppout,'Refine',refine,... 'RelTol',tol,'Abstol',atol); case 'ode23tb' solh = @ode23tb; opt = odeset('OutputFcn',@ppout,'Refine',refine,... 'RelTol',tol,'Abstol',atol); end exist(dfcn); dfh = str2func(dfcn); cflag = 0; [tp,Xp] = feval(solh,dfh,int,w,opt); % [tp,Xp] = ppdp45(dfcn,[0,-1e6],w,ppdisp); set(stopbutt,'enable','off'); aud = get(ppdispa,'UserData'); newhand(1) = aud.line; X1=[z';Xp]; x1 = [X1,[NaN;tp]]; set(newhand(1),'xdata',x1(:,1),'ydata',x1(:,2),'zdata',x1(:,3)); w = z - offs(1)*VV; nstr = get(dud.notice,'string'); nstr(1:4) = nstr(2:5); nstr{5} = 'The second stable trajectory'; set(dud.notice,'string',nstr); set(stopbutt,'enable','on'); [tp,Xp] = feval(solh,dfh,int,w,opt); % [tp,Xp] = ppdp45(dfcn,[0,-1e6],w,ppdisp); set(stopbutt,'enable','off'); aud = get(ppdispa,'UserData'); newhand(2) = aud.line; X2=[z';Xp]; x2 = [X2,[NaN;tp]]; set(newhand(2),'xdata',x2(:,1),'ydata',x2(:,2),'zdata',x2(:,3)); VV=V(:,I(2)); w = z + offs(2)*VV; nstr = get(dud.notice,'string'); nstr(1:4) = nstr(2:5); nstr{5} = 'The first unstable trajectory'; set(dud.notice,'string',nstr); set(stopbutt,'enable','on'); int = [0,1e6]; [tp,Xp] = feval(solh,dfh,int,w,opt); set(stopbutt,'enable','off'); aud = get(ppdispa,'UserData'); newhand(3) = aud.line; X3=[z';Xp]; x3 = [X3,[NaN;tp]]; set(newhand(3),'xdata',x3(:,1),'ydata',x3(:,2),'zdata',x3(:,3)); w = z - offs(2)*VV; nstr = get(dud.notice,'string'); nstr(1:4) = nstr(2:5); nstr{5} = 'The second unstable trajectory'; set(dud.notice,'string',nstr); set(stopbutt,'enable','on'); [tp,Xp] = feval(solh,dfh,int,w,opt); set(stopbutt,'vis','off'); aud = get(ppdispa,'UserData'); newhand(4) = aud.line; X4=[z';Xp]; x4 = [X4,[NaN;tp]]; set(newhand(4),'xdata',x4(:,1),'ydata',x4(:,2),'zdata',x4(:,3)); eqpt = plot(z(1),z(2),... 'color',ecol,... 'markersize',20,... 'marker','.','Erasemode','normal'); dud.solhand = [dud.solhand;newhand(:)]; dud.ephand = [dud.ephand;eqpt]; set(newhand,'color',scol); set(newhand,'erase','normal'); nstr = get(dud.notice,'string'); nstr(1:4) = nstr(2:5); nstr{5} = 'Ready.'; set(dud.notice,'string',nstr); set(gcf,'UserData',dud); elseif strcmp(action,'zoomin') % 'zoomin' is the callback for the Zoomin menu item. It allows the % user to pick a new display rectangle. set(gcf,'WindowButtonDownFcn','pplane81(''zoom'')',... 'WindowButtonUpFcn','1;','inter','on'); set(gca,'inter','on'); dud = get(gcf,'UserData'); nstr = get(dud.notice,'string'); nstr(1:4) = nstr(2:5); nstr{5} = ['Pick a new display rectangle by clicking and ',... 'dragging the mouse, or by clicking on a point.']; set(dud.notice,'string',nstr); elseif strcmp(action,'zoom') disph = gcf; dud = get(disph,'UserData'); axh = dud.axes; aud = get(axh,'UserData'); DY = aud.DY; w = dud.syst.wind; q1 = get(disph,'currentpoint'); p1 = get(axh,'currentpoint'); p1 = p1(1,1:2); rbbox([q1 0 0],q1); p2 = get(axh,'currentpoint'); p2 = p2(1,1:2); if all(abs(p1'-p2')>0.01*DY) a = [p1;p2]; a = [min(a);max(a)]; DY = (a(2,:) - a(1,:))'; else DY = DY/4; a(1) = max(w(1),p1(1)-DY(1)); a(2) = min(w(2),p1(1)+DY(1)); a(3) = max(w(3),p1(2)-DY(2)); a(4) = min(w(4),p1(2)+DY(2)); DY(1) = a(2) - a(1); DY(2) = a(4) - a(3); end WINvect = a(:)'; dud.syst.wind = WINvect; aud.DY = DY; dwind = [WINvect(1); WINvect(3); -WINvect(2); -WINvect(4)]; aud.cwind = dwind - dud.settings.magn*[aud.DY;aud.DY]; set(axh,'UserData',aud); set(disph,'UserData',dud); set(disph,'WindowButtonDownFcn','pplane81(''down'')',... 'WindowButtonUpFcn',''); pplane81('dirfield',disph); ppset = findobj('name','pplane81 Setup'); if isempty(ppset) pplane81('confused'); else sud = get(ppset,'UserData'); sud.c.wind = WINvect; sud.o.wind = WINvect; set(sud.h.wind(1),'string',num2str(WINvect(1))); set(sud.h.wind(2),'string',num2str(WINvect(2))); set(sud.h.wind(3),'string',num2str(WINvect(3))); set(sud.h.wind(4),'string',num2str(WINvect(4))); set(ppset,'UserData',sud); end elseif strcmp(action,'showbar') sbfig = gcbf; domymenu('menubar','toggletoolbar',sbfig); hhset = get(0,'showhiddenhandles'); set(0,'showhiddenhandles','on'); state = get(sbfig,'toolbar'); if strcmp(state,'figure') name = get(sbfig,'name'); fixtb = ['set(gcbo,''state'',''off'');']; set(findobj(sbfig,'tooltipstr','Print'),... 'clickedcallback','pplane81(''print'');'); switch name case 'pplane81 Display' set(findobj(sbfig,'tooltipstr','Zoom Out'),... 'clickedcallback',['pplane81(''zoomback'');' fixtb]); set(findobj(sbfig,'tooltipstr','Zoom In'),... 'clickedcallback',['pplane81(''zoomin'');' fixtb]); delete(findobj(sbfig,'tooltipstr','Rotate 3D')); case 'pplane81 Linearization' delete(findobj(sbfig,'tooltipstr','Zoom Out')); delete(findobj(sbfig,'tooltipstr','Zoom In')); delete(findobj(sbfig,'tooltipstr','Rotate 3D')); otherwise delete(findobj(sbfig,'tooltipstr','Zoom Out')); delete(findobj(sbfig,'tooltipstr','Zoom In')); fud = get(sbfig,'UserData'); switch fud.type case {4, 5} set(findobj(sbfig,'tooltipstr','Rotate 3D'),'vis','on'); otherwise set(findobj(sbfig,'tooltipstr','Rotate 3D'),'vis','off'); end end sbmh = findobj(sbfig,'label','Show &Toolbar'); set(sbmh,'label','Hide &Toolbar','checked','off'); uni = get(0,'units'); ss = get(0,'screensize'); funit = get(sbfig,'units'); set(sbfig,'units',uni); sw = ss(3);sh = ss(4); fpos = get(sbfig,'pos'); if fpos(2)+fpos(4)>sh-40; fpos(2) = sh - fpos(4) -70; set(sbfig,'pos',fpos); end set(sbfig,'units',funit); else sbmh = findobj(sbfig,'label','Hide &Toolbar'); set(sbmh,'label','Show &Toolbar','checked','off'); end set(0,'showhiddenhandles',hhset) elseif strcmp(action,'dall') % 'dall' is the callback for the Erase all graphics objects. disph = gcf; dud = get(disph,'UserData'); notice = dud.notice; kk = find(ishandle(dud.solhand)); delete(dud.solhand(kk)); dud.solhand = []; kk = find(ishandle(dud.ephand)); delete(dud.ephand(kk)); dud.ephand = []; dud.eqpts = []; kk = find(ishandle(dud.contours)); delete(dud.contours(kk)); dud.contours = []; kk = find(ishandle(dud.ics)); delete(dud.ics(kk)); dud.ics = []; if notice ~= 0 set(dud.butt(1),'enable','off'); end set(disph,'UserData',dud); elseif strcmp(action,'dallsol') % 'dallsol' is the callback for the Erase all solutions option. disph = gcf; dud = get(disph,'UserData'); notice = dud.notice; kk = find(ishandle(dud.solhand)); delete(dud.solhand(kk)); dud.solhand = []; if notice ~= 0 set(dud.butt(1),'enable','off'); end set(disph,'UserData',dud); elseif strcmp(action,'dallep') % 'dallep' is the callback for the Erase all equilibrium points option. disph = gcf; dud = get(disph,'UserData'); notice = dud.notice; kk = find(ishandle(dud.ephand)); delete(dud.ephand(kk)); dud.ephand = []; dud.eqpts = []; if notice ~= 0 set(dud.butt(1),'enable','off'); end set(disph,'UserData',dud); elseif strcmp(action,'dalllev') % 'dalllev' is the callback for the Erase all level curves option. disph = gcf; dud = get(disph,'UserData'); notice = dud.notice; kk = find(ishandle(dud.contours)); delete(dud.contours(kk)); dud.contours = []; if notice ~= 0 set(dud.butt(1),'enable','off'); end set(disph,'UserData',dud); elseif strcmp(action,'dallics') % 'dallics' is the callback for the Erase all marked initial cond. disph = gcf; dud = get(disph,'UserData'); notice = dud.notice; kk = find(ishandle(dud.ics)); delete(dud.ics(kk)); dud.ics = []; if notice ~= 0 set(dud.butt(1),'enable','off'); end set(disph,'UserData',dud); elseif strcmp(action,'newton') % Newton's method to find equilibrium points. Iterlimit = 50; Iter=0; % The increment for calculating approximate derivatives. h=.000001; dfcn = input2; zInit = input1; zNext=zInit(:); functionf = zeros(2,1); functionf=feval(dfcn,0,zNext); functionf=functionf(:); % Allow for large/small solutions. errorlim = norm(functionf,inf)*0.000001; while ( (norm(functionf,inf) > errorlim) & (Iter < Iterlimit) ) Iter = Iter + 1; % Now we calculate the jacobian. for jjw=1:2 sav = zNext(jjw); zNext(jjw) = zNext(jjw) + h; functionfhj = feval(dfcn,0,zNext); functionfhj = functionfhj(:); Jacobian(:,jjw) = (functionfhj-functionf)/h; zNext(jjw) = sav; end zNext = zNext - Jacobian\(functionf); functionf = feval(dfcn,0,zNext); functionf=functionf(:); end if Iter > Iterlimit - 1 fLag=[0;0]; else fLag=[1;1]; for j=1:2 sav = zNext(j); zNext(j) = zNext(j) + h; functionfhj = feval(dfcn,0,zNext); functionfhj = functionfhj(:); Jacobian(:,j) = (functionfhj-functionf)/h; zNext(j) = sav; end end output = [zNext,Jacobian,fLag]; elseif strcmp(action,'paraeval') % Replace a parameter with its value in string form. para = deblank(input1); value = input2; value = ['(',value,')']; str = input3; ll = length(str); lp = length(para); lv = length(value); if strcmp(para,str) str = value elseif (ll >= lp+1) k = findstr(para,str); lk = length(k); lopstr = '(+-*/^'; ropstr = ')+-*/^'; s = []; pos = 1; for jj = 1:lk if (((k(jj) == 1)|(find(lopstr == str(k(jj)-1))))... &((k(jj)+lp-1 == ll)|(find(ropstr == str(k(jj) + lp))))) s = [s,str(pos:(k(jj)-1)),value]; pos = k(jj)+lp; end end str = [s,str(pos:ll)]; end output = str; elseif strcmp(action,'settings') % 'settings' is the call back for the Settings menu option. It sets % up the pplane81 Settings window, which allows the user to interactively % change several parameters that govern the behaviour of the program. dud = get(gcf,'UserData'); data.settings = dud.settings; solver = dud.settings.solver; left = 2; nudge = 3; varl = 60; setfig = findobj('name','pplane81 Settings'); if ~isempty(setfig); delete(setfig); end setfig = figure('name','pplane81 Settings',... 'NumberTitle','off',... 'tag','pplane81','vis','off'); pplane81('figdefault',setfig); set(setfig,'menubar','none'); setfr = uicontrol('style','frame'); speedfr = uicontrol('style','frame'); cwfr = uicontrol('style','frame'); ss = uicontrol('style','text','horiz','center',... 'string',['Settings for ',dud.settings.solver]); nstepcall =[ 'data = get(gcf,''UserData'');',... 'ss = max(round(str2num(get(data.nstep,''string''))),0);',... 'data.settings.refine = ss;',... 'set(data.nstep,''string'',num2str(ss));',... 'set(gcf,''UserData'',data);']; nstep = uicontrol('style','text','horiz','left',... 'string','Number of plot steps per computation step: '); data.nstep = uicontrol('style','edit','call',nstepcall,... 'string',num2str(data.settings.refine)); if strcmp(solver,'Runge-Kutta 4') rtolstr = 'Step size: '; rtolentry = num2str(data.settings.stepsize); rtolcall =['data = get(gcf,''UserData'');',... 'data.settings.stepsize = str2num(get(data.rtol,''string''));',... 'set(data.rtol,''string'',num2str(data.settings.stepsize));',... 'set(gcf,''UserData'',data);']; else rtolstr = 'Relative error tolerance: '; rtolentry = num2str(data.settings.tol); rtolcall =['data = get(gcf,''UserData'');',... 'data.settings.tol = str2num(get(data.rtol,''string''));',... 'set(data.rtol,''string'',num2str(data.settings.tol));',... 'set(gcf,''UserData'',data);']; end rtol = uicontrol('style','text','horiz','left',... 'string',rtolstr); data.rtol = uicontrol('style','edit','call',rtolcall,... 'string',rtolentry); kk = 1+2*data.settings.magn; magcall = ['data = get(gcf,''UserData'');',... 'mag = (str2num(get(data.mag,''string''))-1)/2;',... 'data.settings.magn = max(0,mag);',... 'set(gcf,''UserData'',data);']; speedcall = ['data = get(gcf,''UserData'');',... 'me = data.speed;',... 'val = round(get(me,''value''));',... 'set(me,''value'',val);',... 'set(data.sp.val,''string'',num2str(val));',... 'data.settings.speed = val;',... 'set(gcf,''UserData'',data);']; data.speed = uicontrol('style','slider',... 'string','Steps per second.',... 'min',2,... 'max',100,... 'value',data.settings.speed,... 'sliderstep',[1/98,10/98],... 'call',speedcall); data.sp.min = uicontrol('style','text','string',' 2',... 'horiz','left'); data.sp.max = uicontrol('style','text','string','100 ',... 'horiz','right'); data.sp.val = uicontrol('style','text',... 'string',num2str(data.settings.speed),... 'horiz','center'); pps = uicontrol('style','text',... 'string','Solution steps per second.'); cw1 = uicontrol('style','text','horiz','left',... 'string','The calculation window is'); data.mag = uicontrol('style','edit','call',magcall,... 'string',num2str(kk)); cw2 = uicontrol('style','text','horiz','left',... 'string',' times as large as the '); cw3 = uicontrol('style','text','horiz','left',... 'string','display window.'); gob = uicontrol('style','push','string','Go Away','call','close'); chb = uicontrol('style','push',... 'string','Change settings',... 'call','pplane81(''setchange'')'); frsep = 1; ext = get(nstep,'extent'); ht = ext(4)+nudge; stw = ext(3); % = nstepw = rtolw varl = varl*ht/19; cwfrb = 2*left + ht; cwfrh = 2*nudge + 2*ht; speedfrb = cwfrb + cwfrh + frsep; speedfrh = 4*nudge + 3*ht; setfrb = speedfrb + speedfrh + frsep; setfrh = 3*nudge + 3*ht; figh = setfrb + setfrh + left; bb = left; cw3b = cwfrb + nudge; cw2b = cw3b + ht; % = cw1b = cweb rtolb = setfrb + nudge; rtolbb = rtolb; if strcmp(solver,'Runge-Kutta 4') rtolbb = rtolbb +ht/2; end nstepb = rtolb + ht; ssb = nstepb + ht + nudge; ssw = stw + varl; spb1 = speedfrb+2*nudge; spb2 = spb1+ht; spb3 = spb2+ht+nudge; sptw = varl; sl = left + nudge; sptl1 = sl; sptsep = (ssw-3*sptw-sl)/2; sptl2 = sptl1 + sptw + sptsep; sptl3 = sptl2 + sptw + sptsep; cw1ext = get(cw1,'extent'); cw1w = cw1ext(3); cw2ext = get(cw2,'extent'); cw2w = cw2ext(3); cwew = 40*ht/19; frw = 2*nudge + ssw; figw = frw + 2*left; buttw = figw/3; sep = figw/9; sl = left + nudge; gobl = sep; chbl = 2*sep + buttw; sunit = get(0,'units'); set(0,'units','pix'); ssize = get(0,'screensize'); figb = ssize(4) - figh - 50; set(setfig,'pos',[20,figb,figw,figh]); set(setfr,'pos',[left,setfrb,frw,setfrh]); set(speedfr,'pos',[left,speedfrb,frw,speedfrh]); set(cwfr,'pos',[left,cwfrb,frw,cwfrh]); set(ss,'pos',[sl,ssb,ssw,ht]); set(nstep,'pos',[sl,nstepb,stw,ht]); set(data.nstep,'pos',[sl+stw,nstepb,varl,ht]); set(rtol,'pos',[sl,rtolbb,stw,ht]); set(data.rtol,'pos',[sl+stw,rtolbb,varl,ht]); set(data.speed,'pos',[sl,spb1,ssw,ht],... 'backgroundcolor',0.6*[1 1 1],... 'foregroundcolor','r'); set(data.sp.min,'pos',[sptl1,spb2,sptw,ht]); set(data.sp.max,'pos',[sptl3,spb2,sptw,ht]); set(data.sp.val,'pos',[sptl2,spb2,sptw,ht]); set(pps,'pos',[sl,spb3,ssw,ht]); set(cw1,'pos',[sl,cw2b,cw1w,ht]); set(cw2,'pos',[sl+cw1w+cwew,cw2b,cw2w,ht]); set(cw3,'pos',[sl,cw3b,ssw,ht]); set(data.mag,'pos',[sl+cw1w,cw2b,cwew,ht]); set(gob,'pos',[gobl,bb,buttw,ht]); set(chb,'pos',[chbl,bb,buttw,ht]); set(setfig,'UserData',data); set(setfig,'units','normal'); set(findobj(setfig,'type','uicontrol'),'units','normal'); set(setfig,'vis','on','resize','on'); ch = get(setfig,'child'); set([nstep,data.nstep],'vis','off') if strcmp(solver,'Runge-Kutta 4') ch = setdiff(ch,[nstep,data.nstep]); end edith = findobj(setfig,'style','edit'); set(edith,'backgroundcolor','w'); set(ch,'vis','on'); elseif strcmp(action,'setchange') % 'setchange' is the callback for the Change button % on the pplane81 Settings window. data = get(gcf,'UserData'); settings = data.settings; ppdisp = findobj('name','pplane81 Display'); dud = get(ppdisp,'UserData'); if isempty(ppdisp) pplane81('confused'); else dud.settings = settings; set(ppdisp,'UserData',dud); WINvect = dud.wmat; WINvect = WINvect(size(WINvect,1),:); dwind = [WINvect(1); WINvect(3); -WINvect(2); -WINvect(4)]; aud = get(dud.axes,'UserData'); aud.cwind = dwind - dud.settings.magn*[aud.DY;aud.DY]; set(dud.axes,'UserData',aud); end ppset = findobj('name','pplane81 Setup'); if isempty(ppset) pplane81('confused'); else sud = get(ppset,'UserData'); sud.settings = settings; set(ppset,'UserData',sud); end close elseif strcmp(action,'delete') % 'delete' is the callback for the Delete a graphics object selection % on the menu. disph = gcf; dud = get(disph,'UserData'); arr = dud.arr; lv = get(arr.lines,'vis'); av = get(arr.arrows,'vis'); if ~isempty(arr.hx) nv = get(arr.hx(1),'vis'); elseif ~isempty(arr.hy) nv = get(arr.hx(1),'vis'); else nv = zeros(1,0); end if ~isempty(arr.barrows) bv = get(arr.barrows,'vis'); else bv = zeros(1,0); end handles = [arr.lines;arr.arrows;arr.hx;arr.hy;arr.barrows]; set(handles,'vis','off'); oldcall = get(disph,'WindowButtonDownFcn'); set(disph,'WindowButtonDownFcn',''); trjh = dud.solhand; notice = dud.notice; if notice ~= 0 % Display window nstr = get(notice,'string'); nstr(1:4) = nstr(2:5); nstr{5} = 'Select a graphics object with the mouse.'; set(notice,'string',nstr); end % ppginput(1); ginput(1); objh = get(disph,'currentobject'); typ = get(objh,'type'); axh = dud.axes; hh = get(dud.title.axes,'children'); hh = [hh;get(axh,'title');get(axh,'xlabel');get(axh,'ylabel')]; if notice ~= 0 % Display window eph = dud.ephand; levh = dud.contours; eqpt = dud.eqpts; ics = dud.ics; nstr(1:4) = nstr(2:5); if strcmp(typ,'line') eqk = find(eph == objh); if ~isempty(eqk); % An equilibrium point. if ~isempty(eqk) if size(eqpt,1) == 1; dud.eqpts = []; else dud.eqpts = [eqpt(1:eqk-1,:);eqpt(eqk+1:size(eqpt,1),:)]; end end end dud.ics = setdiff(dud.ics,objh); dud.solhand = setdiff(dud.solhand,objh); dud.ephand = setdiff(dud.ephand,objh); dud.contours = setdiff(dud.contours,objh); delete(objh); nstr{5} = 'Ready.'; elseif strcmp(typ,'text') & ~ismember(objh,hh) dud.contours = setdiff(dud.contours,objh)'; delete(objh); nstr{5} = 'Ready.'; else nstr{5} = 'The object you selected cannot be deleted.'; end set(notice,'string',nstr); else % Linearization window if strcmp(typ,'line') eqk = find(dud.ephand == objh); if isempty(eqk); % Npt an equilibrium point. dud.solhand = setdiff(dud.solhand,objh)'; dud.ics = setdiff(dud.ics,objh); delete(objh); end elseif ~ismember(objh,hh) delete(objh); end end set(arr.lines,'vis',lv); set(arr.arrows,'vis',av); set([arr.hx;arr.hy],'vis',nv); set(arr.barrows,'vis',bv); set(disph,'UserData',dud); set(disph,'WindowButtonDownFcn','pplane81(''down'')'); elseif strcmp(action,'text') ppdisp = gcf; oldcall = get(ppdisp,'WindowButtonDownFcn'); set(ppdisp,'WindowButtonDownFcn',''); prompt = ['Enter the text here. Then choose ',... 'the location in the Display Window.']; txtstr = inputdlg(prompt,'Text entry'); if ~isempty(txtstr) txtstr = txtstr{1}; figure(ppdisp); gtext(txtstr); end set(ppdisp,'WindowButtonDownFcn',oldcall); elseif strcmp(action,'plotxy') % Start a graph type = input1; ppdisp =gcf; dud = get(ppdisp,'UserData'); handles = dud.solhand; printstr = dud.printstr; exsol=0; lll = length(handles); if (lll) % Are there any orbits? kk=0; while (kk < lll) & (exsol == 0) kk = kk + 1; xd = get(handles(kk),'xdata'); if (length(xd) > 7) exsol = 1; end end end if (exsol == 0) nstr = get(dud.notice,'string'); nstr(1:3) = nstr(3:5); nstr{4} = 'There are no solution curves.'; nstr{5} = 'Ready.'; set(dud.notice,'string',nstr); return end oldcall = get(ppdisp,'WindowButtonDownFcn'); set(ppdisp,'WindowButtonDownFcn',''); arr = dud.arr; lv = get(arr.lines,'vis'); av = get(arr.arrows,'vis'); if ~isempty(arr.hx) nv = get(arr.hx(1),'vis'); elseif ~isempty(arr.hy) nv = get(arr.hy(1),'vis'); else nv = zeros(1,0); end if ~isempty(arr.barrows) bv = get(arr.barrows,'vis'); else bv = zeros(1,0); end handles = [arr.lines;arr.arrows;arr.hx;arr.hy;arr.barrows]; set(handles,'vis','off'); nstr = get(dud.notice,'string'); nstr(1:4) = nstr(2:5); nstr{5} = 'Select a solution curve with the mouse.'; set(dud.notice,'string',nstr); % ppginput(1); ginput(1); objh = get(ppdisp,'currentobject'); if isempty(objh) solflag = 0; elseif ~(strcmp(get(objh,'type'),'line')) solflag = 0; else t=get(objh,'zdata'); if length(t)<10 solflag = 0; else solflag = 1; end end if solflag == 0 nstr = get(dud.notice,'string'); nstr(1:4) = nstr(2:5); nstr{5} = 'This object is not a solution curve.'; set(dud.notice,'string',nstr); set(ppdisp,'WindowButtonDownFcn',oldcall); set(arr.lines,'vis',lv); set(arr.arrows,'vis',av); set([arr.hx;arr.hy],'vis',nv); return end set(arr.lines,'vis',lv); set(arr.arrows,'vis',av); set([arr.hx;arr.hy],'vis',nv); set(arr.barrows,'vis',bv); t=get(objh,'zdata'); x=get(objh,'xdata'); y=get(objh,'ydata'); hh = pplane81('plotxyfig',type,gcf); set(ppdisp,'WindowButtonDownFcn',oldcall); aud = get(hh(2),'UserData'); ud = get(hh(1),'UserData'); ud.data = [t;x;y]; set(hh(1),'UserData',ud); pplane81('plxy',aud.rad,hh(1)); elseif strcmp(action,'plotxyfig') % Build the graph figure. type = input1; cfig = input2; % This is ppdisp or a ppxy ppdisp = findobj('name','pplane81 Display'); dud = get(ppdisp,'UserData'); ppdispa = dud.axes; xstring = get(get(ppdispa,'XLabel'),'string'); ystring = get(get(ppdispa,'YLabel'),'string'); dud = get(ppdisp,'UserData'); ppxy = figure('number','off',... 'tag','pplane81','visible','off'); name = ['pplane81 t-plot # ',int2str(ppxy)]; ud.color = dud.color; ud.type = 0; set(ppxy,'name',name,'UserData',ud); pplane81('figdefault',ppxy); % Menus and Toolbar hhsetup = get(0,'showhiddenhandles'); set(0,'showhiddenhandles','on'); % Configure the Toolbar. set(ppxy,'ToolBar','none'); % Menus emenu = findobj(gcf,'label','&Edit'); tmenu = findobj(gcf,'label','&Tools'); delete(tmenu) % File menu fmenu = findobj(gcf,'label','&File'); delete(findobj(fmenu,'label','&New Figure')); delete(findobj(fmenu,'label','&Open...')); set(findobj(fmenu,'label','&Close'),'pos',1); set(findobj(fmenu,'label','&Save'),... 'pos',2,'separator','off'); set(findobj(fmenu,'label','Save &As...'),... 'pos',3); set(findobj(fmenu,'label','&Export...'),... 'pos',4); delete(findobj(fmenu,'label','Pre&ferences...')); set(findobj(fmenu,'label','Pa&ge Setup...'),'pos',4); set(findobj(fmenu,'label','Print Set&up...'),'pos',5); set(findobj(fmenu,'label','Print Pre&view...'),'pos',6); set(findobj(fmenu,'label','&Print...'),'pos',7); qmenu = uimenu(fmenu,'label','Quit pplane81',... 'call','pplane81(''quit'')',... 'separator','on'); % Insert menu imenu = findobj(gcf,'label','&Insert'); inschild = get(imenu,'child'); legitem = findobj(inschild,'label','&Legend'); colitem = findobj(inschild,'label','&Colorbar'); delete([legitem,colitem]); % View menu set(findobj(gcf,'label','&Figure Toolbar'),... 'call','pplane81(''showbar'')'); set(0,'showhiddenhandles',hhsetup); fs = get(ppxy,'defaultaxesfontsize'); r = fs/10; axw = 437*.8; % Default axes width axh = 315*.8; % Default axes height axl = 45; % Default axes left buttw = 40; % Default button width titleh = 45; % Default title height legw = 90; % Default legend width axb = 35; % Default axes bottom sep = 10; fh = axb + axh + titleh; % Default figure height. fw = axl + axw + sep + legw + sep; % Default figure width. uni = get(0,'units'); set(0,'units','pixels'); ss = get(0,'screensize'); set(0,'units',uni); sw = ss(3);sh = ss(4); if r*fh > sh -80; r = (sh-80)/fh; end if r*fw > sw - 50 r = (sw - 50)/fw; end fs = 10*r; lw = 0.5*r; set(ppxy,'defaultaxesfontsize',fs,'defaultaxeslinewidth',lw); set(ppxy,'defaulttextfontsize',fs); set(ppxy,'defaultlinelinewidth',lw); set(ppxy,'defaultuicontrolfontsize',fs*0.9); axw = r*axw; axh = r*axh; axl = r*axl; legw = r*legw; axb = r*axb; sep = r*sep; % The legend. leg = axes('units','pix',... 'box','on',... 'drawmode','fast',... 'nextplot','add',... 'xtick',[-1],'ytick',[-1],... 'xticklabel','','yticklabel','',... 'xlim',[0,1],'ylim',[0,1],... 'clipping','on','visible','off'); set(leg,'UserData',str2mat(xstring,ystring)); axes(leg); xh = text(0,0,xstring,'units','norm','visible','off',... 'parent',leg); yh = text(0,0,ystring,'units','norm','visible','off',... 'parent',leg); set(xh,'units','pix') set(yh,'units','pix') xext = get(xh,'extent'); yext = get(yh,'extent'); set(xh,'units','norm') set(yh,'units','norm') th = max(xext(4),yext(4))+3*r; axhn = max(axh, 11*th + 4 + 6*r); % Nudge the axes height if needed. axw = axw*axhn/axh; axh = axhn; frsep = (axh - 11*th - 4 - 6*r)/4; legh = 2*(th+3*r); legb = axb + axh - legh; %frsep*4 + th*9 +4; legl = axl + axw + sep; % If the legend text is too big, make the legend and the figure larger. legpos = [legl+1, legb, legw, legh]; set(leg,'pos',legpos); xext = get(xh,'extent'); yext = get(yh,'extent'); tw = max(xext(3),yext(3))+0.1; rrrr = 0.5; if (tw > rrrr) rrr = tw/rrrr; legw = legw*rrr; tw = rrrr; end legpos = [legl+1, legb, legw, legh]; set(leg,'pos',legpos); fw = axl + axw + sep + legw + sep; % Figure width. tx = min(1-tw,2/3); set(xh,'pos',[tx,2/3]); set(yh,'pos',[tx,1/3]); line('xdata',[0.1,tx-0.1],'ydata',[2/3,2/3],'linestyle','-',... 'color',ud.color.tx,'vis','off'); line('xdata',[0.1,tx-0.1],'ydata',[1/3,1/3],'linestyle','--',... 'color',ud.color.ty,'vis','off'); % Set up the title. tax = axes('box','off',... 'xlim',[0 1],'ylim',[0 1],... 'units','pix','vis','off',... 'xtick',[-1],'ytick',[-1],... 'xticklabel','','yticklabel',''); titb = axb + axh; titl = axl; titw = axw; titeq = text(0.01,0.5,dud.tstr,'vert','middle'); set(titeq,'units','pix'); ext = get(titeq ,'extent'); set(titeq,'units','data'); titleh = ext(4)+6*r; p1 = get(dud.title.p1,'string'); p2 = get(dud.title.p2,'string'); p3 = get(dud.title.p3,'string'); pp1 = text(0.85,0.5,p1,'vert','middle'); ext = get(pp1,'extent'); pp1l = 1 - ext(3); set(pp1,'pos',[pp1l,0.5]); pp2 = text(0.7,0.5,p2,'vert','middle'); ext = get(pp2,'extent'); pp2l = pp1l - ext(3)-0.05; set(pp2,'pos',[pp2l,0.5]); pp3 = text(0.7,0.5,p3,'vert','middle'); ext = get(pp3,'extent'); pp3l = pp2l - ext(3)-0.05; set(pp3,'pos',[pp3l,0.5]); fh = axb + axh + titleh; % Figure height. pos = get(cfig,'pos'); fl = max(0,pos(1)-30); fb = max(20,pos(2)-30); fpos = [fl,fb,fw,fh]; set(ppxy,'pos',fpos); set(tax,'pos',[titl,titb,titw,titleh]); set(tax,'color',get(ppxy,'color')); % The t-plot axes. axpos = [axl,axb,axw,axh]; axy = axes('units','pix','pos',axpos,'box',... 'on','xgrid','on','ygrid','on','next','add',... 'drawmode','fast'); output = [ppxy,axy]; % The position of the close button bbot = axb; bwid = legw; bleft = legl; bpos = [bleft,bbot,bwid,th]; % The position of the print button. pbbot = bbot + th + frsep; pbpos = [bleft, pbbot, legw, th]; % The position of the crop button. cbbot = pbbot + th + frsep; cbpos = [bleft, cbbot, legw, th]; % The positions of the radio frame and its elements. frbot = cbbot + th + frsep; frleft = legl; frht = 6*th+4; frpos = [frleft,frbot,legw,frht]; frtitleft = frleft+2; frtitwid = legw-4; frtitpos = [frtitleft, frbot+5*th, frtitwid, th]; radleft = frleft + 2; radwid = legw - 4; radbot5 = frbot + 2; radbot4 = radbot5 + th; radbot3 = radbot4 + th; radbot2 = radbot3 + th; radbot1 = radbot2 + th; pcall = [ 'ppdisp = findobj(''name'',''pplane81 Display'');',... 'dud = get(ppdisp,''UserData'');',... 'eval(dud.printstr)' ]; but = uicontrol('style','push',... 'pos',bpos,... 'string','Go away',... 'call','close'); pbut = uicontrol('style','push',... 'pos',pbpos,... 'string','Print',... 'call',pcall); cbut = uicontrol('style','push',... 'pos',cbpos,... 'string','Crop',... 'call','pplane81(''crop'')'); pframe = uicontrol('style','frame','pos',frpos); pfrtitle = uicontrol('style','text','pos',frtitpos,... 'string','Graph','horizon','center'); rval1 = ~(type -1); rval2 = (~(type -2))*2; rval3 = (~(type -3))*3; rval4 = (~(type -4))*4; rval5 = (~(type -5))*5; rad(1) = uicontrol('style','radio',... 'pos',[radleft radbot1 radwid th],... 'string',xstring,'value',rval1); rad(2) = uicontrol('style','radio',... 'pos',[radleft radbot2 radwid th],... 'string',ystring,'value',rval2,'max',2); rad(3) = uicontrol('style','radio',... 'pos',[radleft radbot3 radwid th],... 'string','Both','value',rval3,'max',3); rad(4) = uicontrol('style','radio',... 'pos',[radleft radbot4 radwid th],... 'string','3 D','value',rval4,'max',4); rad(5) = uicontrol('style','radio',... 'pos',[radleft radbot5 radwid th],... 'string','Composite','value',rval5,'max',5); for i=1:5 set(rad(i),'UserData',[rad(:,[1:(i-1),(i+1):5]),leg,axy]); end callrad = [ 'me = gcbo;',... 'vv = get(me,''value'');'... 'mm = get(me,''max'');'... 'if vv ,'... ' hand = get(me,''UserData'');',... ' set(hand(1:4),''value'',0);',... ' pplane81(''plxy'',me,gcf);'... 'end, '... 'set(me,''value'',mm);'... 'axy = gca;',... 'aud = get(axy,''UserData'');',... 'aud.rad = me;',... 'set(axy,''UserData'',aud);']; set(rad,'call',callrad); set(findobj(ppxy,'type','axes'),'units','normal'); set(findobj(ppxy,'type','uicontrol'),'units','normal'); set(ppxy,'resize','on'); set(ppxy,'visible','on'); nstr = get(dud.notice,'string'); nstr(1:4) = nstr(2:5); nstr{5} = 'Ready.'; set(dud.notice,'string',nstr); aud.h = []; aud.int = [0 0]; aud.crop = cbut; aud.rad = rad(type); aud.leg = leg; set(axy,'UserData',aud); set(cbut,'UserData',axy); elseif strcmp(action,'plxy') radbut = input1; fig = input2; fud = get(fig,'UserData'); figure(fig) tbh = findobj(allchild(fig),'flat','type','uitoolbar'); r3dh = findobj(tbh,'tooltipstr','Rotate 3D'); ud = get(radbut,'UserData'); axy = ud(6); axis('auto'); delete(get(axy,'children')); leg = ud(5); legch = get(leg,'children'); type = get(radbut,'max'); fud.type = type; set(fig,'UserData',fud); if type == 3 set(leg,'visible','on'); set(legch,'visible','on'); else set(leg,'visible','off'); set(legch,'visible','off'); end aud = get(axy,'UserData'); aud.h = []; set(axy,'UserData',aud); strings = get(leg,'UserData'); xstring = deblank(strings(1,:)); ystring = deblank(strings(2,:)); ud = get(fig,'UserData'); data = ud.data; color = ud.color; t = data(1,:); x = data(2,:); y = data(3,:); tmin = min(t); tmax = max(t); xmin = min(x); xmax = max(x); ymin = min(y); ymax = max(y); et = max((tmax-tmin)/20,1e-4); ex = max((xmax-xmin)/20,1e-4); ey = max((ymax-ymin)/20,1e-4); tmin1 = tmin - et; tmax1 = tmax + et; xmin = xmin - ex; xmax = xmax + ex; ymin = ymin - ey; ymax = ymax + ey; axes(axy); switch type case 1 plot(t,x,'color',color.tx,'linestyle','-'); view(2); axis([tmin,tmax,xmin,xmax]) xlabel('t') ylabel(xstring) set(r3dh,'vis','off') set(fig,'WindowButtonDownFcn','pplane81(''plxybdf'')'); case 2 plot(t,y,'color',color.ty,'linestyle','-'); view(2); axis([tmin,tmax,ymin,ymax]) xlabel('t') ylabel(ystring) set(r3dh,'vis','off') rotate3d off set(fig,'WindowButtonDownFcn','pplane81(''plxybdf'')'); case 3 plot(t,x,'color',color.tx,'linestyle','-'); plot(t,y,'color',color.ty,'linestyle','--'); axis([tmin,tmax,min(xmin,ymin),max(xmax,ymax)]) view(2); xlabel('t') ylabel([xstring,' and ',ystring]) set(r3dh,'vis','off') rotate3d off set(fig,'WindowButtonDownFcn','pplane81(''plxybdf'')'); case 4 plot3(x,y,t,'color',color.tx,'linestyle','-'); view([-30,20]); axis([xmin,xmax,ymin,ymax,tmin1,tmax1]); zlabel('t') xlabel(xstring); ylabel(ystring); set(r3dh,'vis','on') rotate3d off set(fig,'WindowButtonDownFcn',''); otherwise ax = [xmin,xmax,ymin,ymax,tmin1,tmax1]; view([-30,20]); axis(ax); plot3(x,y,ax(5)*ones(size(x)),'color',color.ty,'linestyle','-'); plot3(x,ax(4)*ones(size(y)),t,'color',color.ty,'linestyle','-'); plot3(ax(2)*ones(size(x)),y,t,'color',color.ty,'linestyle','-'); plot3(x,y,t,'color',color.tx,'linestyle','-'); zlabel('t') xlabel(xstring); ylabel(ystring); set(r3dh,'vis','on') rotate3d off set(fig,'WindowButtonDownFcn',''); end if (type < 4) set(aud.crop,'vis','on','enable','off'); else set(aud.crop,'vis','off'); end if type < 4 set(fig,'windowbuttondownfcn','pplane81(''plxybdf'')'); elseif type == 5 set(fig,'windowbuttondownfcn',' '); end elseif strcmp(action,'crop') cb = gcbo; axy = get(cb,'UserData'); aud = get(axy,'UserData'); delete(aud.h); set(aud.crop,'enable','off'); aud.h = []; set(axy,'UserData',aud); ppxy = gcf; ud = get(ppxy,'UserData'); data = ud.data; color = ud.color; t = data(1,:); x = data(2,:); y = data(3,:); int = aud.int; k = find((t>=int(1)) & (t<=int(2))); t = t(k); x = x(k); y = y(k); type = get(aud.rad,'max'); hh = pplane81('plotxyfig',type,ppxy); ud = get(hh(1),'UserData'); ud.data = [t;x;y]; set(hh(1),'UserData',ud); aud = get(hh(2),'UserData'); pplane81('plxy',aud.rad,hh(1)) elseif strcmp(action,'plxybdf') dispa = gca; dispf = gcf; set(dispf,'windowbuttonmotionfcn','pplane81(''plxybmf'')',... 'windowbuttonupfcn','pplane81(''plxybuf'')',... 'inter','on'); aud = get(dispa,'UserData'); delete(aud.h); aud.h = []; set(aud.crop,'enable','off'); point = get(dispa,'currentpoint'); aud.start = point(1,1); aud.finish = point(1,1); % aud.centh = plot(point(1,1),point(1,2),'or','erase','xor'); aud.h = plot(point(1,1),point(1,2),'--g',... 'erase','xor','vis','off'); set(dispa,'UserData',aud); elseif strcmp(action,'plxybmf') dispa = gca; aud = get(dispa,'UserData'); point = get(dispa,'currentpoint'); finish = point(1,1); start = aud.start; xlim = get(dispa,'xlim'); ylim = get(dispa,'ylim'); if abs(finish - start)>0.05*(xlim(2)-xlim(1)); set(aud.h,'xdata',[start,start,NaN,finish,finish],... 'ydata',[ylim,NaN,ylim],... 'vis','on'); end aud.finish = finish; set(gca,'UserData',aud); elseif strcmp(action,'plxybuf') dispa = gca; dispf = gcf; xlim = get(dispa,'xlim'); ylim = get(dispa,'ylim'); aud = get(dispa,'UserData'); set(gcf,'windowbuttonmotionfcn','',... 'windowbuttonupfcn','',... 'inter','on'); start = aud.start; finish = aud.finish; if abs(finish - start)>0.05*(xlim(2)-xlim(1)); set(aud.h,'erase','normal','vis','on'); set(aud.crop,'enable','on'); aud.int = [min(start,finish),max(start,finish)]; else set(aud.h,'vis','off'); end set(gca,'UserData',aud); elseif strcmp(action,'print') dud = get(gcf,'UserData'); nstr = get(dud.notice,'string'); nstr(1:4) = nstr(2:5); nstr{5} = 'Preparing to print the pplane81 Display Window. Please be patient.'; set(dud.notice,'string',nstr); nstr(1:4) = nstr(2:5); nstr{5} = 'Printing the pplane81 Display Window.'; set(dud.notice,'string',nstr); set(gcf,'pointer','watch'); eval(dud.printstr); nstr(1:4) = nstr(2:5); nstr{5} = 'Ready.'; set(dud.notice,'string',nstr); set(gcf,'pointer','arrow'); elseif strcmp(action, 'eqptlist') % The labels for equilibrium points. disp(' '); EqPtType = ['Saddle point. '; 'Nodal sink. '; 'Nodal source. '; 'Spiral sink. '; 'Spiral source. '; 'Spiral equilibrium point.'; 'Source. '; 'Sink. '; 'Unspecified. ';]; dud = get(gcf,'UserData'); EqPtList = dud.eqpts; if isempty(EqPtList) disp('No equilibrium points have been computed.'),disp(' ') else disp('The following equilibrium points have been calculated:') disp(' ') L = size(EqPtList,1); eqpttext = cell(L,1); for k = 1:L disp([sprintf('(%6.4f, %6.4f)\t',EqPtList(k,1),EqPtList(k,2)),... EqPtType(EqPtList(k,3),:)]) end disp(' ') end elseif strcmp (action,'zoomback') disph = gcf; dud = get(disph,'UserData'); axh = dud.axes; Xname = dud.syst.xvar; Yname = dud.syst.yvar; wmat = dud.wmat; WINvect = dud.syst.wind; NN = size(wmat,1); wch = 0;j=0; while wch == 0 j = j+1; if WINvect == wmat(j,:) wch = j; end end winstr = cell(1,NN); for j = 1:NN a = num2str(wmat(j,1)); b = num2str(wmat(j,2)); c = num2str(wmat(j,3)); d = num2str(wmat(j,4)); winstr{j} = [' ',a,' < ',Xname,' < ',b,' & ',c,' < ', Yname,' < ',d]; end [sel,ok] = listdlg('liststring',winstr,... 'selectionmode','single',... 'listsize',[400,200],... 'initialvalue',wch,... 'name','pplane81 Zoomback',... 'promptstring','Select a rectangle:',... 'OKString','Zoom'); if (~isempty(sel)) WINvect = wmat(sel,:); dud.syst.wind = WINvect; set(gcf,'UserData',dud); ppset = findobj('name','pplane81 Setup'); sud = get(ppset,'UserData'); set(disph,'UserData',dud); set(sud.h.wind(1),'string',num2str(WINvect(1))); set(sud.h.wind(2),'string',num2str(WINvect(2))); set(sud.h.wind(3),'string',num2str(WINvect(3))); set(sud.h.wind(4),'string',num2str(WINvect(4))); sud.c.wind = WINvect; sud.o.wind = WINvect; set(ppset,'UserData',sud); aud = get(axh,'UserData'); aud.DY = [WINvect(2) - WINvect(1);WINvect(4) - WINvect(3)]; dwind = [WINvect(1); WINvect(3); -WINvect(2); -WINvect(4)]; aud.cwind = dwind - dud.settings.magn*[aud.DY;aud.DY]; set(axh,'UserData',aud); pplane81('dirfield',disph); end elseif strcmp(action,'figdefault') fig = input1; set(fig,'CloseRequestFcn','pplane81(''closefcn'')'); ppset = findobj('name','pplane81 Setup'); sud = get(ppset,'UserData'); ud = get(fig,'UserData'); ud.ssize = sud.ssize; fs = sud.fontsize; ud.fontsize = fs; style = sud.style; set(fig,'defaulttextfontsize',fs); set(fig,'defaultaxesfontsize',fs); set(fig,'defaultuicontrolfontsize',9*fs/10) lw = 0.5*fs/10; set(fig,'defaultaxeslinewidth',lw) set(fig,'defaultlinelinewidth',lw) set(fig,'defaultaxesfontname','helvetica') set(fig,'defaultaxesfontweight','normal') switch style case 'black' % if isunix | isvms, gamma = 0.5; else gamma = 0.0; end whitebg(fig,[0,0,0]) if isunix | isvms fc = [.35 .35 .35]; else fc = [.2 .2 .2]; end set(fig,'color',fc); set(fig,'defaultaxescolor',[0 0 0]) % whitebg(fig,brighten([.2 .2 .2],gamma)) set(fig,'defaultaxescolor',[0 0 0]) set(fig,'defaultaxescolororder', ... 1-[0 0 1;0 1 0;1 0 0;0 1 1;1 0 1;1 1 0;.25 .25 .25]) % ymcrgbw % set(fig,'colormap',brighten(jet(64),gamma)) set(fig,'colormap',jet(64)) set(fig,'defaultsurfaceedgecolor',[0 0 0]); case 'white' whitebg(fig,[1 1 1]) set(fig,'color',[.8 .8 .8]) set(fig,'defaultaxescolor',[1 1 1]) set(fig,'defaultaxescolororder', ... [0 0 1;0 .5 0;1 0 0;0 .75 .75;.75 0 .75;.75 .75 0;.25 .25 .25]) % bgrymck set(fig,'colormap',jet(64)) set(fig,'defaultsurfaceedgecolor',[0 0 0]) case 'display' whitebg(fig,[1 1 1]) set(fig,'defaultaxescolor',[1 1 1]) set(fig,'defaultaxescolororder', ... [0 0 1;0 .5 0;1 0 0;0 .75 .75;.75 0 .75;.75 .75 0;.25 .25 .25]) % bgrymck set(fig,'colormap',jet(64)) set(fig,'defaultsurfaceedgecolor',[0 0 0]) set(fig,'color',[1 1 1]*240/255); set(fig,'defaultuicontrolbackgroundcolor',[1 1 1]*220/255); set(fig,'defaultaxesfontweight','bold') set(fig,'defaulttextfontweight','bold') set(fig,'defaultaxeslinewidth',1) set(fig,'defaultlinelinewidth',1) case 'bw' whitebg(fig,[0 0 0]) set(fig,'color',[0 0 0]) set(fig,'defaultaxescolor',[0 0 0]) set(fig,'defaultaxescolororder', ... [1 1 1]) set(fig,'colormap',[1 1 1;0 0 0]) set(fig,'defaultsurfaceedgecolor',[0 0 0]) end set(fig,'UserData',ud); elseif strcmp(action,'zoominsq') % 'zoominsq' is the callback for the Zoom in square menu item. set(gcf,'WindowButtonDownFcn','pplane81(''zoomsqd'')',... 'WindowButtonUpFcn','1;','inter','on'); set(gca,'inter','on'); dud = get(gcf,'UserData'); nstr = get(dud.notice,'string'); nstr(1:4) = nstr(2:5); nstr{5} = ['Pick a center and ',... 'drag the mouse, or just click on a center.']; set(dud.notice,'string',nstr); elseif strcmp(action,'zoomsqd') ppdispa = gca; aud = get(ppdispa,'UserData'); point = get(ppdispa,'currentpoint'); aud.center = point(1,[1,2]); aud.centh = plot(point(1,1),point(1,2),'or','erase','xor'); aud.box = plot(point(1,1),point(1,2),'--r','erase','xor'); set(ppdispa,'UserData',aud); set(gcf,'windowbuttonmotionfcn','pplane81(''zoomsqm'')',... 'windowbuttonupfcn','pplane81(''zoomsqu'')'); elseif strcmp(action,'zoomsqm') ppdispa = gca; aud = get(ppdispa,'UserData'); point = get(ppdispa,'currentpoint'); point = point(1,[1,2])'; un = get(ppdispa,'units'); set(ppdispa,'units','pix'); w = get(ppdispa,'pos'); w = w([3,4]);w = w(:); set(ppdispa,'units',un); cent = aud.center(:); lam = max(abs(point - cent)./w); v = lam*w; data = cent(:,[1 1 1 1 1]) + [v,[-v(1);v(2)], -v, [v(1);-v(2)],v]; set(aud.box,'xdata',data(1,:),'ydata',data(2,:)); elseif strcmp(action,'zoomsqu') disph = gcf; dud = get(disph,'UserData'); axh = dud.axes; aud = get(axh,'UserData'); hand = [aud.centh;aud.box]; set(hand,'erase','normal'); delete(hand); DY = aud.DY(:); un = get(axh,'units'); set(axh,'units','pix'); w = dud.syst.wind; point = get(axh,'currentpoint'); point = point(1,[1,2]);point = point(:); cent = aud.center(:); ww = get(axh,'pos'); ww = ww([3,4]); ww = ww(:); lamb = max(abs(point-cent)./ww); v = lamb*ww; if ~all(v > 0.01*DY) points = [w([1 1 2 2]);w([3 4 3 4])]; for j=1:4 lambs(j) = max(abs(points(:,j)-cent)./ww); end lamb = min(lambs); v = lamb*ww/4; end p1 = cent + v; p2 = cent - v; a = [p2';p1']; a = [min(a);max(a)]; DY = (a(2,:) - a(1,:))'; WINvect = a(:)'; dud.syst.wind = WINvect; aud.DY = DY; dwind = [WINvect(1); WINvect(3); -WINvect(2); -WINvect(4)]; aud.cwind = dwind - dud.settings.magn*[aud.DY;aud.DY]; set(axh,'units',un,'UserData',aud); set(disph,'UserData',dud); set(disph,'WindowButtonDownFcn','pplane81(''down'')',... 'WindowButtonMotionFcn','pplane81(''cdisp'')',... 'WindowButtonUpFcn',''); pplane81('dirfield',disph); ppset = findobj('name','pplane81 Setup'); if isempty(ppset) pplane81('confused'); else sud = get(ppset,'UserData'); sud.c.wind = WINvect; sud.o.wind = WINvect; set(sud.h.wind(1),'string',num2str(WINvect(1))); set(sud.h.wind(2),'string',num2str(WINvect(2))); set(sud.h.wind(3),'string',num2str(WINvect(3))); set(sud.h.wind(4),'string',num2str(WINvect(4))); set(ppset,'UserData',sud); end elseif strcmp(action,'level') ppdisp = gcf; dud = get(ppdisp,'UserData'); lfcn = dud.level; Xname = dud.syst.xvar; Yname = dud.syst.yvar; Xname(find(abs(Xname)==92))=[]; % Remove \s if any. Yname(find(abs(Yname)==92))=[]; pplevel = findobj('name','pplane81 Level sets'); if ~isempty(pplevel) delete(pplevel) end pplevel = figure('name','pplane81 Level sets',... 'vis','off',... 'NumberTitle','off','tag','pplane81'); pplane81('figdefault',pplevel); set(pplevel,'menubar','none'); lev.fr1 = uicontrol('style','frame'); lev.fr2 = uicontrol('style','frame'); inst1str = ['Enter the function in terms of the variables ',... Xname, ' and ', Yname,':']; lev.inst1 = uicontrol('style','text','horiz','left',... 'string',inst1str); lev.lfcn = uicontrol('style','edit','horiz','center',... 'string',lfcn,'call','',... 'background','w'); lev.ch(3) = uicontrol('style','radio','horiz','center',... 'min',0,'max',3,... 'value',0,... 'vis','on',... 'string','Let pplane81 decide.'); lev.ch(2) = uicontrol('style','radio','horiz','center',... 'min',0,'max',2,... 'value',0,... 'string','Select a point in the Display Window.'); lev.inst2 = uicontrol('style','text','horiz','left',... 'string',['Choose one of the following ways to',... ' choose level value(s):']); lev.ch(1) = uicontrol('style','radio','horiz','center',... 'min',0,'max',1,... 'value',0,... 'string','Enter a vector of level values.'); lev.rhs = uicontrol('style','edit','horiz','center',... 'string',' ','call',''); lev.proc = uicontrol('style','push',... 'string','Proceed',... 'call','pplane81(''levcomp'')'); lev.close = uicontrol('style','push',... 'string','Close',... 'call','close'); for i=1:3 set(lev.ch(i),'UserData',lev.ch(:,[1:(i-1),(i+1):3])); end callrad = [ 'me = get(gcf,''currentobject'');',... 'kk = get(me,''max'');',... 'col = get(me,''backg'');',... 'set(get(me,''UserData''),''value'',0),',... 'set(me,''value'',kk);',... 'ud = get(gcf,''UserData'');',... 'if kk == 1,',... ' set(ud.rhs,''enable'',''on'',''backg'',''w'');',... 'else,',... ' set(ud.rhs,''enable'',''off'',''backg'',col);',... 'end,']; set(lev.ch,'call',callrad); left = 2; varl = 300; buttw = 60; nudge = 3; tab = 15; lines1 = 5; lines2 = 2; xex = get(lev.inst1,'extent'); ht = xex(4)+nudge; frw = varl + 2*tab+ 2*nudge; fr1bot = 2*left + ht; fr1ht = lines1*(nudge + ht) + nudge; fr2bot = fr1bot + fr1ht; fr2ht = lines2*(nudge + ht) + nudge; vbot = fr1bot + nudge; ch1bot = vbot + nudge + ht; ch2bot = ch1bot + nudge + ht; ch3bot = ch2bot + nudge + ht; inst2bot = ch3bot + nudge + ht; fbot = fr2bot + nudge; inst1bot = fbot + nudge + ht; fleft = left + nudge + tab; instleft = left + nudge; chleft = instleft + tab; vleft = chleft + tab; vw = (frw - 2*vleft); fr1wind = [left,fr1bot,frw,fr1ht]; fr2wind = [left,fr2bot,frw,fr2ht]; inst1wind = [instleft,inst1bot,varl,ht]; inst2wind = [instleft,inst2bot,varl,ht]; fwind = [fleft,fbot,varl,ht]; ch1wind = [chleft,ch1bot,varl,ht]; ch2wind = [chleft,ch2bot,varl,ht]; ch3wind = [chleft,ch3bot,varl,ht]; vwind = [vleft,vbot,vw,ht]; figw = 2*left + frw; fight = 3*left + ht + fr1ht + fr2ht; figwind = [40, 300, figw, fight]; buttw = frw/2; sep = (figw - 2*buttw)/3; closel = sep; procl = 2*sep+buttw; clwind = [closel,left,buttw,ht]; procwind = [procl,left,buttw,ht]; set(pplevel,'pos',figwind); set(lev.fr1,'pos',fr1wind); set(lev.fr2,'pos',fr2wind); set(lev.inst1,'pos',inst1wind); set(lev.inst2,'pos',inst2wind); set(lev.ch(1),'pos',ch1wind); set(lev.ch(2),'pos',ch2wind); set(lev.ch(3),'pos',ch3wind); set(lev.rhs,'pos',vwind); set(lev.lfcn,'pos',fwind); set(lev.proc,'pos',procwind); set(lev.close,'pos',clwind); set(lev.ch(3),'value',3); set(lev.rhs,'enable','off'); child = get(pplevel,'children'); set(pplevel,'vis','on','UserData',lev); set(child,'vis','on'); elseif strcmp(action,'levcomp') pplevel = gcf; ud = get(pplevel,'UserData'); ppdisp = findobj('name','pplane81 Display'); dud = get(ppdisp,'UserData'); ppset = findobj('name','pplane81 Setup'); sud = get(ppset,'UserData'); ch = ud.ch; val = zeros(1,3); for kk = 1:3 val(kk) = get(ch(kk),'value'); end KK = max(val); lfcn = get(ud.lfcn,'string'); l=length(lfcn); for ( k = fliplr(findstr('.',lfcn))) if (find('*/^' == lfcn(k+1))) lfcn = [lfcn(1:k-1), lfcn(k+1:l)]; end l=l-1; end pnameh = sud.h.pname; pvalh = sud.h.pval; pflag = zeros(1,4); perr = []; lfcn(find(abs(lfcn)==32))=[]; for kk = 1:4; pn = get(pnameh(kk),'string'); pv = get(pvalh(kk),'string'); if ~isempty(pn) pn(find(abs(pn)==92))=[]; if isempty(pv) perr = pvalh(kk); else pv(find(abs(pv)==32))=[]; lfcn = pplane81('paraeval',pn,pv,lfcn); end end end l = length(lfcn); for (k=fliplr(find((lfcn=='^')|(lfcn=='*')|(lfcn=='/')))) lfcn = [lfcn(1:k-1) '.' lfcn(k:l)]; l = l+1; end WINvect = dud.syst.wind; XxXxXx = WINvect(1) + rand*(WINvect(2)-WINvect(1)); YyYyYy = WINvect(3) + rand*(WINvect(4)-WINvect(3)); Xname = dud.syst.xvar; Yname = dud.syst.yvar; Xname(find(abs(Xname)==92))=[]; % Remove \s if any. Yname(find(abs(Yname)==92))=[]; err = 0;res = 1; eval([Xname,'=XxXxXx;'],'err = 1;'); eval([Yname,'=YyYyYy;'],'err = 1;'); eval(['res = ',lfcn, ';'],'err = 1;'); if err | isempty(res) errmsg = 'The function does not evaluate correctly.'; fprintf('\a') errordlg(errmsg,'PPLANE error','on'); return; end Xmin = WINvect(1); Xmax = WINvect(2); Ymin = WINvect(3); Ymax = WINvect(4); N = 50; k = 4; deltax=(Xmax - Xmin)/(N-1); deltay=(Ymax - Ymin)/(N-1); XXXg=Xmin + deltax*[-k:N+k]; YYYg=Ymin + deltay*[-k:N+k]; [Xx,Yy]=meshgrid(XXXg,YYYg); Xxx=Xx(:);Yyy=Yy(:); Ww = zeros(size(Xxx)); eval([Xname,'=Xxx;'],'err = 1;'); eval([Yname,'=Yyy;'],'err = 1;'); eval(['Ww = ',lfcn, ';']) KKK = 3; %# of significant figures. switch KK case 1 % vector input rhs = get(ud.rhs,'string'); rhs = str2num(rhs); case 2 % mouse input figure(ppdisp); % [XX,YY] = ppginput(1); [XX,YY] = ginput(1); figure(pplevel); eval([Xname,'=XX;'],'err = 1;'); eval([Yname,'=YY;'],'err = 1;'); eval(['rhs = ',lfcn, ';'],'err = 1'); LL = ceil(log10(abs(rhs))); rhs = round(10^(KKK-LL)*rhs); rhs = 10^(LL-KKK)*rhs; case 3 % pplane81 input Www = Ww; kkk = find(isnan(Www)); Www(kkk) = []; kkk = find(imag(Www)); Www(kkk) = []; MM = max(Www); mm = min(Www); LL = ceil(log10(MM-mm)); NN = 7; % Number of curves rhs = mm+(1:NN).^2*(MM-mm)/NN^2; rhs = round(10^(KKK-LL)*rhs); rhs = 10^(LL-KKK)*rhs; end Ww = reshape(Ww,N+2*k+1,N+2*k+1); lrhs = length(rhs); if lrhs == 0 return elseif lrhs == 1 rhs = [rhs,rhs]; end figure(ppdisp); [Cm,hcont] = contour(Xx,Yy,Ww,rhs,'--'); hlabel = clabel(Cm,hcont); % set(hlabel,'fontsize',dud.fontsize,... % 'color',dud.color.level,... % 'rotation',0); set(hlabel,'fontsize',dud.fontsize,... 'color',[1,0,0],... 'rotation',0); set(hcont,'visible','on',... 'color',dud.color.level,... 'linestyle',':'); dud.contours = [dud.contours ;hcont;hlabel]; set(ppdisp,'UserData',dud); elseif strcmp(action,'restart') ppset = findobj('name','pplane81 Setup'); oldfiles = dir('pptp*.m'); for k = 1:length(oldfiles) fn = oldfiles(k).name; fid = fopen(fn,'r'); ll = fgetl(fid); ll = fgetl(fid); ll = fgetl(fid); fclose(fid); if strcmp(ll,'%% Created by pplane81') delete(fn) end end h = findobj('tag','pplane81'); delete(setdiff(h,ppset)); sud = get(ppset,'UserData'); sud.flag = 0; set(ppset,'UserData',sud); figure(ppset) elseif strcmp(action,'quit') ppset = findobj('name','pplane81 Setup'); sud = get(ppset,'UserData'); if sud.remtd rmpath(tempdir); end oldfiles = dir([tempdir,'pptp*.m']); for k = 1:length(oldfiles) fn = [tempdir,oldfiles(k).name]; fid = fopen(fn,'r'); ll = fgetl(fid); ll = fgetl(fid); ll = fgetl(fid); fclose(fid); if strcmp(ll,'%% Created by pplane81') delete(fn) end end h = findobj('tag','pplane81'); delete(h); elseif strcmp(action,'closefcn') fig = gcf; name = get(fig,'name'); if strcmp(name,'pplane81 Setup') | strcmp(name,'pplane81 Display') quest = ['Closing this window will cause all pplane81 ',... 'windows to close, and pplane81 will stop. ',... 'Do you want to quit pplane81?']; butt = questdlg(quest,'Quit pplane81?','Quit','Cancel','Quit'); if strcmp(butt,'Quit') pplane81('quit'); end elseif strcmp(name,'pplane81 Linearization') dud = get(fig,'UserData'); fcn = dud.function; if (exist(fcn)==2) delete([fcn,'.m']);end delete(findobj('label',name)); delete(fig); else delete(findobj('label',name)); delete(fig); end elseif strcmp(action,'confused') tstring = 'pplane81 is totally confused'; qstring = {['You will have to restart pplane81 from '... 'the beginning in order to ',... 'do anything new. However, it might be possible '... 'to save the current system ',... 'or the gallery to make your restart easier, '... 'or it may be possible to ',... 'print out a figure, if the appropriate '... 'figures are visible. In such a case ',... 'it would be best to do nothing now.']; 'What do you want to do?'}; bstr1 = 'Quit and restart pplane81.'; bstr2 = 'Just quit pplane81.'; bstr3 = 'Do nothing.'; answer = questdlg(qstring,tstring,bstr1,bstr2,bstr3,bstr1); if strcmp(answer,bstr1) delete(findobj('tag','pplane81')); pplane81;return elseif strcmp(answer,bstr2) delete(findobj('tag','pplane81')); return else return end elseif strcmp(action,'export') % export is the callback for the Export solution data item in the % Options menu. disph = gcf; dud = get(disph,'UserData'); arr = dud.arr; lv = get(arr.lines,'vis'); av = get(arr.arrows,'vis'); if ~isempty(arr.hx) nv = get(arr.hx(1),'vis'); elseif ~isempty(arr.hy) nv = get(arr.hx(1),'vis'); else nv = zeros(1,0); end if ~isempty(arr.barrows) bv = get(arr.barrows,'vis'); else bv = zeros(1,0); end handles = [arr.lines;arr.arrows;arr.hx;arr.hy;arr.barrows]; set(handles,'vis','off'); oldcall = get(disph,'WindowButtonDownFcn'); set(disph,'WindowButtonDownFcn',''); trjh = dud.solhand; notice = dud.notice; switch length(trjh) case 0 if notice ~= 0 nstr = get(notice,'string'); nstr(1:3) = nstr(3:5); nstr{4} = 'There are no solutions.'; nstr{5} = 'Ready.'; set(notice,'string',nstr); end th = []; case 1 th = trjh; otherwise if notice ~= 0 nstr = get(notice,'string'); nstr(1:4) = nstr(2:5); nstr{5} = 'Select a solution with the mouse.'; set(notice,'string',nstr); end % ppginput(1); ginput(1); th = get(disph,'currentobject'); end if isempty(th) if notice ~= 0 nstr = get(notice,'string'); nstr(1:3) = nstr(3:5); nstr{4} = 'The item selected is not a solution.'; nstr{5} = 'Ready.'; set(notice,'string',nstr); end else vars = evalin('base','who'); no = 1; kk = 0; while no kk = kk + 1; vstr = ['ppdata',num2str(kk)]; if ~any(strcmp(vars,vstr)) no = 0; end end yname = dud.syst.yvar; if abs(yname(1)) == 92 yname = yname(2:length(yname)); end xname = dud.syst.xvar; if abs(xname(1)) == 92 xname = xname(2:length(xname)); end tname = 't'; tval = get(th,'zdata'); tval = tval(:); xval = get(th,'xdata'); xval = xval(:); yval = get(th,'ydata'); yval = yval(:); ivstr = struct(tname,tval,xname,xval,yname,yval); assignin('base',vstr,ivstr); if notice ~= 0 nstr = get(notice,'string'); nstr(1:3) = nstr(3:5); nstr{4} = ['The data has been exported as the structure ',... vstr,' with fields ', tname, ', ', xname, ', and ', yname,'.']; nstr{5} = 'Ready.'; set(notice,'string',nstr); end end set(arr.lines,'vis',lv); set(arr.arrows,'vis',av); set([arr.hx;arr.hy],'vis',nv); set(arr.barrows,'vis',bv); set(disph,'UserData',dud); set(disph,'WindowButtonDownFcn','pplane81(''down'')'); elseif strcmp(action,'cdisp') [ppcbo,ppdisp] = gcbo; dud = get(ppdisp,'UserData'); cp = get(ppdisp,'currentpoint'); fpos = get(ppdisp,'pos'); ppax = dud.axes; xd = get(ppax,'xlim'); yd = get(ppax,'ylim'); apos = get(ppax,'pos'); xp = xd(1) + (cp(1) - apos(1)*fpos(3))*(xd(2)-xd(1))/(apos(3)*fpos(3)); yp = yd(1) + (cp(2) - apos(2)*fpos(4))*(yd(2)-yd(1))/(apos(4)*fpos(4)); str = ['(',num2str(xp,3),', ',num2str(yp,3),')']; set(dud.ccwind,'string',str); elseif strcmp(action,'periodic') % Find a periodic orbit. ppdisp = findobj('name','pplane81 Display'); dud = get(ppdisp,'UserData'); dfcn = dud.function; direction = dud.dir; notice = dud.notice; settings = dud.settings; refine = settings.refine; tol = settings.tol; AA = -1e6; BB = 1e6; switch direction case 0 intplus = [0, BB]; intminus = [0, AA]; case -1 intplus = [0, 0]; intminus = [0, AA]; case 1 intplus = [0, BB]; intminus = [0, 0]; end ppdispa = dud.axes; aud = get(ppdispa,'UserData'); DY = aud.DY; aud.plot = 0; atol = tol*DY*1e-4'; set(ppdispa,'UserData',aud); nstr = get(dud.notice,'string'); nstr(1:4) = nstr(2:5); nstr{5} = 'Choose a starting point with the mouse.'; set(dud.notice,'string',nstr); % z00 = ppginput(1); z00 = ginput(1); z0 = z00; ptstr = [' (',num2str(z0(1),2), ', ', num2str(z0(2),2), ')']; solver = dud.settings.solver; opt = odeset('OutputFcn',@ppout,'Refine',refine,... 'RelTol',tol,'Abstol',atol); switch solver case 'Dormand Prince' solh = @ppdp45; opt = ppdisp; case 'Runge-Kutta 4' solh = @pprk4; opt = ppdisp; case 'ode45' solh = @ode45; case 'ode23' solh = @ode23; case 'ode113' solh = @ode113; case 'ode15s' solh = @ode15s; case 'ode23s' solh = @ode23s; case 'ode23t' solh = @ode23t; case 'ode23tb' solh = @ode23tb; end exist(dfcn); dfh = str2func(dfcn); dud.noticeflag = 0; % Notices only from here. set(ppdisp,'UserData',dud); if intplus(2)>intplus(1) if notice ~= 0 nstr = get(notice,'string'); nstr(1:4) = nstr(2:5); nstr{5} = 'The forward orbit '; set(notice,'string',nstr); end drawnow [tp,xp] = feval(solh,dfh,intplus,z0,opt); aud = get(ppdispa,'UserData'); set(aud.line,'erase','normal') delete(aud.line); aud.line = []; z0 = aud.y; stop = aud.stop; nstr = get(dud.notice,'string'); switch stop case 1 nstr{5} = [nstr{5}, ' left the computation window.']; case 2 zz = aud.zz; ystr = ['(',num2str(zz(1),2), ', ', num2str(zz(2),2),').']; nstr{5} = [nstr{5}, ' --> a possible eq. pt. near ',ystr]; case 3 nstr{5} = [nstr{5}, ' --> a closed orbit']; case 4 nstr{5} = [nstr{5}, ' was stopped by the user.']; case 5 ystr = ['(',num2str(y(1),2), ', ', num2str(y(2),2),').']; nstr(1:3) = nstr(2:4); nstr{4} = [nstr{5},' experienced a failure at ',ystr]; nstr{5} = 'Problem is singular or tolerances are too large.'; end set(dud.notice,'string',nstr); drawnow if stop==3 % periodic orbit found. % [tp,xp] = feval(solh,dfh,intplus,z0,opt); rr = sqrt(((xp(:,1) - z0(1))/DY(1)).^2 + ((xp(:,2) - z0(2))/DY(2)).^2); kk = find(rr < 0.01); LL = length(kk); kmax = max(kk); kkk = max(kk(find(kk0.0001*abs(T) & kkk < 10 [tt,yy] = feval(solh,dfh,[0,T],z0,opt); yy = yy(length(tt),:); dyy = feval(dfh,0,yy'); NN = ((yy-z0')*dz0)/(dyy'*dz0); T = T - NN; kkk = kkk +1; end nstr{5} = [nstr{5}, ' of period ', num2str(abs(T),3), '.']; set(dud.notice,'string',nstr); drawnow [tp,xp] = feval(solh,dfh,[0 3*T],z0,opt); aud.gstop = 1; set(ppdispa,'UserData',aud); hnew = plot(xp(:,1),xp(:,2),'color','b'); set(hnew,'zdata',tp); solhand = [dud.solhand;hnew]; dud.solhand = solhand; end drawnow end if intminus(2) a possible eq. pt. near ',ystr]; case 3 nstr{5} = [nstr{5}, ' --> a nearly closed orbit.']; case 4 nstr{5} = [nstr{5}, ' was stopped by the user.']; case 5 ystr = ['(',num2str(y(1),2), ', ', num2str(y(2),2),').']; nstr(1:3) = nstr(2:4); nstr{4} = [nstr{5},' experienced a failure at ',ystr]; nstr{5} = 'Problem is singular or tolerances are too large.'; end set(dud.notice,'string',nstr); drawnow if stop==3 % periodic orbit found. z0 = aud.y; rr = sqrt(((xp(:,1) - z0(1))/DY(1)).^2 + ((xp(:,2) - z0(2))/DY(2)).^2); kk = find(rr < 0.01); LL = length(kk); kmax = max(kk); kkk = max(kk(find(kk0.0001*abs(T) & kkk < 10 [tt,yy] = feval(solh,dfh,[0,T],z0,opt); yy = yy(length(tt),:); dyy = feval(dfh,0,yy'); NN = ((yy-z0')*dz0)/(dyy'*dz0); T = T - NN; kkk = kkk +1; end nstr{5} = [nstr{5}, ' of period ', num2str(abs(T),3), '.']; set(dud.notice,'string',nstr); drawnow [tp,xp] = feval(solh,dfh,[0 3*T],z0,opt); aud.gstop = 1; set(ppdispa,'UserData',aud); hnew = plot(xp(:,1),xp(:,2),'color','b'); set(hnew,'zdata',tp); solhand = [dud.solhand;hnew]; dud.solhand = solhand; end drawnow end dud.noticeflag = 1; dud.dir = 0; set(ppdisp,'UserData',dud); aud.plot = 1; set(ppdispa,'UserData',aud); if notice ~= 0 nstr = get(notice,'string'); nstr(1:4) = nstr(2:5); nstr{5} = 'Ready'; set(notice,'string',nstr); end drawnow end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [tout,yout] = ppdp45(dfcn,tspan,y0,disph) % PPDP45 is an implementation of the explicit Runge-Kutta (4,5) which % is described in Chapters 5 & 6 of John Dormand's book, % Numerical Methods for Differential Equations. % % This is the same algorithm used in ODE45, part of the new MATLAB % ODE Suite. Details are to be found in The MATLAB ODE Suite, % L. F. Shampine and M. W. Reichelt, SIAM Journal on Scientific % Computing, 18-1, 1997. % Input the user data. dud = get(disph,'UserData'); dispha = dud.axes; ud = get(dispha,'UserData'); refine = dud.settings.refine; tol = dud.settings.tol; gstop = ud.gstop; plotf = ud.plot; DY = ud.DY; col = dud.color.temp; speed = dud.settings.speed; slow = (speed < 100); % Initialize the stopping criteria. if gstop % Initialize detection of closed orbits & limit cycles. % Choose a random direction & initialize the search for orbit maxima in % that direction. theta = rand*2*pi; R = [cos(theta), sin(theta); -sin(theta),cos(theta)]; qq = R*(y0(:)); rr = [qq,qq]; z = DY(1) + sqrt(-1)*DY(2); w=exp(i*theta); a1 = w*z; a2 = w*(z'); a=max(abs(real([a1,a2]))); b=max(abs(imag([a1,a2]))); perpeps = a*0.00001; % a/2000; % We want to be real % close in this direction. paraeps = b/100; % We do not have to be % so careful in this direction. tk = 0; turn = zeros(2,10); % The test for an equilibrium point. sinkeps = 0.005/refine; % The compute window. cwind = ud.cwind; end % The stop button. stop = 0; ud.stop = 0; % Set the the line handle. ph = plot([y0(1),y0(1)],[y0(2),y0(2)],... 'color',col,... 'erase','none',... 'parent',dispha); ud.line = ph; set(dispha,'UserData',ud); % Initialize the loop. t0 = tspan(1); tfinal = tspan(2); tdir = sign(tfinal - t0); t = t0; y = y0(:); % By default, hmax is 1/10 of the interval. hmax = min(abs(0.1*(tfinal-t)),1); rDY = DY(:,ones(1,refine)); steps = 0; block = 120; tout = zeros(block,1); yout = zeros(block,2); N = 1; tout(N) = t; yout(N,:) = y.'; % Initialize method parameters. pow = 1/5; % C = [1/5; 3/10; 4/5; 8/9; 1; 1]; % Not needed because the sytem is autonomous. A = [ 1/5 3/40 44/45 19372/6561 9017/3168 0 9/40 -56/15 -25360/2187 -355/33 0 0 32/9 64448/6561 46732/5247 0 0 0 -212/729 49/176 0 0 0 0 -5103/18656 0 0 0 0 0 0 0 0 0 0 ]; bhat = [35/384 0 500/1113 125/192 -2187/6784 11/84 0]'; % E = bhat - b. E = [71/57600; 0; -71/16695; 71/1920; -17253/339200; 22/525; -1/40]; if refine > 1 sigma = (1:refine-1) / refine; S = cumprod(sigma(ones(4,1),:)); bstar = [ 1 -183/64 37/12 -145/128 0 0 0 0 0 1500/371 -1000/159 1000/371 0 -125/32 125/12 -375/64 0 9477/3392 -729/106 25515/6784 0 -11/7 11/3 -55/28 0 3/2 -4 5/2 ]; bstar = bstar*S; end f = zeros(2,7); f0 = feval(dfcn,t,y); mm = max(abs(f0./DY)); absh = hmax; if mm absh = min(absh,1/(100*mm)); end f(:,1) = f0; minNsteps =20; % THE MAIN LOOP tic while ~stop % hmin is a small number such that t+h is distinquishably % different from t if abs(h) > hmin. hmin = 16*eps*abs(t); absh = min(hmax, max(hmin, absh)); h = tdir * absh; % LOOP FOR ADVANCING ONE STEP. while stop~=5 % hC= h * C; hA = h * A; f(:,2) = feval(dfcn, t, y + f*hA(:,1)); f(:,3) = feval(dfcn, t, y + f*hA(:,2)); f(:,4) = feval(dfcn, t, y + f*hA(:,3)); f(:,5) = feval(dfcn, t, y + f*hA(:,4)); f(:,6) = feval(dfcn, t, y + f*hA(:,5)); tn = t + h; yn = y + f*h*bhat; f(:,7) = feval(dfcn, tn, yn); % Estimate the error. err = abs(h * f * E); alpha = (2*max(err./((abs(y)+abs(yn)+DY*1e-3)*tol)))^pow; if alpha < 1 % Step is OK break else if absh <= hmin % This keeps us out of an infinite loop. stop = 5; break; end absh = max(hmin,0.8*absh / alpha); h = tdir * absh; end % if alpha < 1 end % while stop~=5 steps = steps + 1; oldN = N; N = N + refine; if N > length(tout) tout = [tout; zeros(block,1)]; yout = [yout; zeros(block,2)]; end if refine > 1 % computed points, with refinement j = oldN+1:N-1; tout(j) = t + h*sigma'; yout(j,:) = (y(:,ones(length(j),1)) + f*h*bstar).'; tout(N) = tn; yout(N,:) = yn.'; else % computed points, no refinement tout(N) = tn; yout(N,:) = yn.'; end % Update stop. Maybe the stop button has been pressed. ud = get(dispha,'UserData'); stop = max(ud.stop,stop); if gstop % Are we out of the compute window? yl = yout(N,:).'; if any([yl;-yl] - cwind < 0); stop = 1; end % If the step in the phase plane is small we assume there is a sink. if (steps > minNsteps) yy = yout(oldN:N,:).'; dyy = yy(:,1:refine) - yy(:,2:(refine+1)); dyy = dyy./rDY; MMf = min(sqrt(sum(dyy.^2))); if (MMf<=sinkeps*absh); z0 = yy(:,refine+1); zz = pplane81('newton',z0,dfcn); zz = zz(:,1); ud.zz = zz; nz = norm((zz-z0)./DY); if nz <= 0.01; stop = 2; ud.y = z0; end minNsteps = minNsteps + 20; end end % We look for a maximum in the randomly chosen direction. If % we find one, we compare it with previously found maxima. If % our new one is close to an old one we stop. If not, we % record the on. jj = oldN+1:N; yy = yout(jj,:).'; rrr = R*yy; v = [rr,rrr]; rr = v(:,[refine+1,refine+2]); % Use this next time. [m,ii] = max(v(1,:)); if( 1< min(ii) & max(ii) size(turn,2) turn = [turn,zeros(2,10)]; end turn(:,tk+1) = v(:,ii); end elseif (abs(tn-tfinal) < hmin) stop = 6; end % if gstop if plotf ii = oldN:N; set(ph,'Xdata',yout(ii,1),'Ydata',yout(ii,2)); drawnow end % if plotf % Compute a new step size. absh = max(hmin,0.8*absh / max( alpha,0.1)); absh = min(absh,tdir*(tfinal - tn)); h = absh*tdir; % Advance the integration one step. t = tn; y = yn; f(:,1) = f(:,7); % Already evaluated % dfcn(tnew,ynew) if slow ttt= N/(speed*refine); tt = toc; while tt < ttt tt = toc; end end end % while ~stop ud.stop = stop; set(dispha,'UserData',ud); tout = tout(1:N); yout = yout(1:N,:); if dud.notice ~= 0 && dud.noticeflag nstr = get(dud.notice,'string'); switch stop case 1 nstr{5} = [nstr{5}, ' left the computation window.']; case 2 ystr = ['(',num2str(zz(1),2), ', ', num2str(zz(2),2),').']; nstr{5} = [nstr{5}, ' --> a possible eq. pt. near ',ystr]; case 3 nstr{5} = [nstr{5}, ' --> a nearly closed orbit.']; case 4 nstr{5} = [nstr{5}, ' was stopped by the user.']; case 5 ystr = ['(',num2str(y(1),2), ', ', num2str(y(2),2),').']; nstr(1:3) = nstr(2:4); nstr{4} = [nstr{5},' experienced a failure at ',ystr]; nstr{5} = 'Problem is singular or tolerances are too large.'; end set(dud.notice,'string',nstr); drawnow end %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [tout,yout] = pprk4(dfcn,tspan,y0,disph) % PPRK4 is an implementation of the fourth order Runge-Kutta method. % Input the user data. dud = get(disph,'UserData'); dispha = dud.axes; ud = get(dispha,'UserData'); refine = dud.settings.refine; tol = dud.settings.tol; gstop = ud.gstop; ssize = dud.settings.stepsize; plotf = ud.plot; DY = ud.DY; col = dud.color.temp; speed = dud.settings.speed; slow = (speed < 100); % Initialize the stopping criteria. % The compute window. if gstop % Initialize detection of closed orbits & limit cycles. % Choose a random direction & initialize the search for orbit maxima in % that direction. theta = rand*2*pi; R = [cos(theta), sin(theta); -sin(theta),cos(theta)]; qq = R*(y0(:)); rr = [qq,qq]; z = DY(1) + sqrt(-1)*DY(2); w=exp(i*theta); a1 = w*z; a2 = w*(z'); a=max(abs(real([a1,a2]))); b=max(abs(imag([a1,a2]))); perpeps = a*0.00001; % We want to be real % close in this direction. paraeps = b/100; % We do not have to be % so careful in this direction. tk = 0; turn = zeros(2,10); % The test for an equilibrium point. sinkeps = 0.0001; % The compute window. cwind = ud.cwind; end % The stop button. stop = 0; ud.stop = 0; % Set the the line handle. ph = plot([y0(1),y0(1)],[y0(2),y0(2)],'color',col,... 'erase','none',... 'parent',dispha); ud.line = ph; set(dispha,'UserData',ud); % Initialize the loop. t0 = tspan(1); tfinal = tspan(2); tdir = sign(tfinal - t0); t = t0; y = y0(:); h = ssize*tdir; steps = 0; block = 120; tout = zeros(block,1); yout = zeros(block,2); N = 1; tout(N) = t; yout(N,:) = y.'; minNsteps =20; % The main loop tic while ~stop if abs(t - tfinal) < ssize h = tfinal - t; end % Compute the slope s1 = feval(dfcn,t,y); s1=s1(:); s2 = feval(dfcn, t + h/2, y + h*s1/2); s2=s2(:); s3 = feval(dfcn, t + h/2, y + h*s2/2); s3=s3(:); s4 = feval(dfcn, t + h, y + h*s3); s4=s4(:); t = t + h; y = y + h*(s1 + 2*s2 + 2*s3 +s4)/6; if N >= length(tout) tout = [tout;zeros(block,1)]; yout = [yout;zeros(block,2)]; end oldN = N; N = N + 1; tout(N) = t; yout(N,:) = y.'; steps = steps + 1; % Update stop. Maybe the stop button has been pressed. ud = get(dispha,'UserData'); stop = max(ud.stop,stop); if gstop % Are we out of the compute window? yl = yout(N,:).'; if any([yl;-yl] - cwind < 0); stop = 1; end % If the step in the phase plane is small we assume there is a sink. if (steps > minNsteps) yy = yout(N-1:N,:).'; dyy = yy(:,1) - yy(:,2); dyy = dyy./DY; MMf = sqrt(sum(dyy.^2)); if (MMf<=sinkeps) z0 = yy(:,refine+1); zz = pplane81('newton',z0,dfcn); zz = zz(:,1); ud.zz = zz; nz = norm((zz-z0)./DY); if nz <= 0.01; stop = 2; ud.y = z0; end minNsteps = minNsteps + 20; end end % We look for a maximum in the randomly chosen direction. If % we find one, we compare it with previously found maxima. If % our new one is close to an old one we stop. If not, we % record the on. yy = yout(N,:).'; rrr = R*yy; v = [rr,rrr]; rr = v(:,[2,3]); % Use this next time. [m,ii] = max(v(1,:)); if( 1< min(ii) & max(ii)<3 ) % If the max is in the middle. kk=0; while ( (kk size(turn,2) turn = [turn,zeros(2,10)]; end turn(:,tk+1) = v(:,ii); end end % if gstop if (abs(t-tfinal) < 0.001*ssize) stop = 6; end if plotf nn = (N-1):N; set(ph,'Xdata',yout(nn,1),'Ydata',yout(nn,2)); drawnow end % if plotf if slow ttt= N/(speed); tt = toc; while tt < ttt tt = toc; end end end % while ~stop ud.stop = stop; set(dispha,'UserData',ud); tout = tout(1:N); yout = yout(1:N,:); if dud.notice ~= 0 && dud.noticeflag nstr = get(dud.notice,'string'); switch stop case 1 nstr{5} = [nstr{5}, ' left the computation window.']; case 2 ystr = ['(',num2str(zz(1),2), ', ', num2str(zz(2),2),').']; nstr{5} = [nstr{5}, ' --> a possible eq. pt. near ',ystr]; case 3 nstr{5} = [nstr{5}, ' --> a nearly closed orbit.']; case 4 nstr{5} = [nstr{5}, ' was stopped by the user.']; case 5 ystr = ['(',num2str(y(1),2), ', ', num2str(y(2),2),').']; nstr(1:3) = nstr(2:4); nstr{4} = [nstr{5},' experienced a failure at ',ystr]; nstr{5} = 'Problem is singular or tolerances are too large.'; end set(dud.notice,'string',nstr); end drawnow %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function output = ppout(t,y,flag,varargin) % PPOUT is the output function for pplane81. % Copywright (c) John C. Polking, Rice University % Last Modified: January 21, 2001 output = 0; ppdisp = findobj(get(0,'child'),'flat','name','pplane81 Display'); %gca; % dud = get(ppdisp,'UserData'); ppdispa = dud.axes; dfcn = dud.function; ud = get(ppdispa,'UserData'); stop = ud.stop; gstop = ud.gstop; col = dud.color.temp; plotf = ud.plot; speed = dud.settings.speed; DY = ud.DY; slow = (speed < 100); if (nargin < 3) | (isempty(flag)) if stop output = 1; vers = version; vers = str2num(vers(1:3)); if vers<6.5 feval(@ppout,t,y,'done'); end else L = length(t); if gstop % Update stop. Are we out of the compute window? yl = y(:,L); if any([yl;-yl] - ud.cwind < 0); stop = 1; end % If the derivative function is small we assume there is a % sink. minNsteps = ud.minNsteps; if (ud.i > minNsteps) yy = [ud.y,y]; dyy = yy(:,1:L) - yy(:,2:(L+1)); rDY = DY(:,ones(1,L)); dyy = dyy./rDY; MMf = min(sqrt(sum(dyy.^2))); if (MMf<=ud.sinkeps*abs(t(1) - t(L))) z0 = yy(:,L+1); zz = pplane81('newton',z0,dfcn); zz = zz(:,1); ud.zz = zz; nz = norm((zz-z0)./DY); if nz <= 0.01; ud.zz = zz; stop = 2; end ud.minNsteps = minNsteps + 20; end else ud.i = ud.i + 1; end % We look for a maximum in the randomly chosen direction. If % we find one, we compare it with previously found maxima. If % our new one is close to an old one we stop. If not, we % record the position. rr = ud.R*y; v = [ud.rr,rr]; ud.rr = v(:,[L+1,L+2]); [m,ii] = max(v(1,:)); %ii = ii(1); if( 1< min(ii) & max(ii)= size(turn,2) ud.turn = [turn,zeros(2,10)]; end ud.turn(:,tk+1) = v(:,ii); end end output = 0; ud.stop = stop; yold = ud.y; ud.y = y(:,L); set(ppdispa,'UserData',ud); if slow ttt = clock; newtime = (24*ttt(4)+ttt(5))*60 + ttt(6); ctime = ud.ctime; N = ud.i; while newtime < ctime + N/speed ttt = clock; newtime = (24*ttt(4)+ttt(5))*60 + ttt(6); end end % Finally we plot the newest line segment. if plotf set(ud.line,'Xdata',[yold(1),y(1,:)],'Ydata',[yold(2),y(2,:)]); drawnow end end else switch(flag) case 'init' % ppout(tspan,y0,'init') if slow ttt = clock; ctime = (24*ttt(4)+ttt(5))*60 + ttt(6); ud.ctime = ctime; end ud.y = y(:); ud.i = 1; % Set the the line handle. figure(ppdisp); ud.line = plot([ud.y(1),ud.y(1)],[ud.y(2),ud.y(2)],... 'color',col,'erase','none'); if gstop % Chose a random direction & initialize the search for orbit % maxima in that direction. theta = rand*2*pi; ud.R = [cos(theta), sin(theta); -sin(theta),cos(theta)]; qq = ud.R*y; ud.rr = [qq,qq]; z = ud.DY(1) + sqrt(-1)*ud.DY(2); w = exp(i*theta); r = abs(z); a1 = w*z; a2 = w*(z'); a = max(abs(real([a1,a2]))); b = max(abs(imag([a1,a2]))); ud.perpeps = a*0.00001;; % We want to be real % close in this direction. ud.paraeps = b/100; % We do not have to be % too careful in this direction. ud.sinkeps = 0.005/dud.settings.refine; ud.minNsteps = 20; ud.tk = 0; ud.turn = zeros(2,10); end output = 0; ud.stop = 0; set(ppdispa,'UserData',ud); case 'done' % ppn6(t,y,'done'); if dud.noticeflag nstr = get(dud.notice,'string'); if ~isempty(y) set(ud.line,'Xdata',[ud.y(1),y(1,:)],'Ydata',[ud.y(2),y(2,:)]); end switch stop case 1 nstr{5} = [nstr{5}, ' left the computation window.']; case 2 yy = ud.zz; ystr = ['(',num2str(yy(1),2), ', ', num2str(yy(2),2),').']; nstr{5} = [nstr{5}, ' --> a possible eq. pt. near ',ystr]; case 3 nstr{5} = [nstr{5}, ' --> a nearly closed orbit.']; case 4 nstr{5} = [nstr{5}, ' was stopped by the user.']; end set(dud.notice,'string',nstr); drawnow output = 1; end end end