clear all;
y=[1.6];
Q=[4.8];
yinv2=1./(y.^2);
B=[2];
dz=[0.0:0.01:2];
ysize=size(y,2);
Qsize=size(Q,2);
Bsize=size(B,2);
dzsize=size(dz,2);
            hold off
            plot(0,0);
            hold on
y2sub=zeros(1,dzsize(1));
y2sup=zeros(1,dzsize(1));
y1sub=zeros(1,dzsize(1));
y1sup=zeros(1,dzsize(1));
E11up=zeros(1,dzsize(1));
E2down=zeros(1,dzsize(1));
for i=1:ysize,
    y1=y(i);
    for j=1:Qsize,
        Q1=Q(j);
        for k=1:Bsize,
            B1=B(k);
            B2=B1;
            Q2=Q1;
            E1=y1+Q1^2/(2*9.81*B1^2*y1^2);
            yc=E1/1.5;
            for l=1:dzsize,
                dz2=dz(l);
                E2=E1-dz(l);
                y2c=(Q2^2/(B2^2*9.81))^(1/3);
                E2c=1.5*y2c;
                dzm=E1-E2c;
                if E2 > E2c
                    E11=E1;
                else
                    E11=E2c+dz2;
                end
                y101=1/600000000000/B1*((-5504587155963302752293577981651332*Q1^2+8000000000000000000000000000000000*E11^3*B1^2+36*Q1*(23379999812960001496319988029439720749282234005742127954062976369*Q1^2-67957866123003737682636765205572000000000000000000000000000000000*E11^3*B1^2)^(1/2))*B1)^(1/3)+200000000000/3*E11^2*B1/((-5504587155963302752293577981651332*Q1^2+8000000000000000000000000000000000*E11^3*B1^2+36*Q1*(23379999812960001496319988029439720749282234005742127954062976369*Q1^2-67957866123003737682636765205572000000000000000000000000000000000*E11^3*B1^2)^(1/2))*B1)^(1/3)+1/3*E11;
                y102=-1/1200000000000/B1*((-5504587155963302752293577981651332*Q1^2+8000000000000000000000000000000000*E11^3*B1^2+36*Q1*(23379999812960001496319988029439720749282234005742127954062976369*Q1^2-67957866123003737682636765205572000000000000000000000000000000000*E11^3*B1^2)^(1/2))*B1)^(1/3)-100000000000/3*E11^2*B1/((-5504587155963302752293577981651332*Q1^2+8000000000000000000000000000000000*E11^3*B1^2+36*Q1*(23379999812960001496319988029439720749282234005742127954062976369*Q1^2-67957866123003737682636765205572000000000000000000000000000000000*E11^3*B1^2)^(1/2))*B1)^(1/3)+1/3*E11+1/200000000000*i*3^(1/2)*(1/6/B1*((-5504587155963302752293577981651332*Q1^2+8000000000000000000000000000000000*E11^3*B1^2+36*Q1*(23379999812960001496319988029439720749282234005742127954062976369*Q1^2-67957866123003737682636765205572000000000000000000000000000000000*E11^3*B1^2)^(1/2))*B1)^(1/3)-20000000000000000000000/3*E11^2*B1/((-5504587155963302752293577981651332*Q1^2+8000000000000000000000000000000000*E11^3*B1^2+36*Q1*(23379999812960001496319988029439720749282234005742127954062976369*Q1^2-67957866123003737682636765205572000000000000000000000000000000000*E11^3*B1^2)^(1/2))*B1)^(1/3));
                y103=-1/1200000000000/B1*((-5504587155963302752293577981651332*Q1^2+8000000000000000000000000000000000*E11^3*B1^2+36*Q1*(23379999812960001496319988029439720749282234005742127954062976369*Q1^2-67957866123003737682636765205572000000000000000000000000000000000*E11^3*B1^2)^(1/2))*B1)^(1/3)-100000000000/3*E11^2*B1/((-5504587155963302752293577981651332*Q1^2+8000000000000000000000000000000000*E11^3*B1^2+36*Q1*(23379999812960001496319988029439720749282234005742127954062976369*Q1^2-67957866123003737682636765205572000000000000000000000000000000000*E11^3*B1^2)^(1/2))*B1)^(1/3)+1/3*E11-1/200000000000*i*3^(1/2)*(1/6/B1*((-5504587155963302752293577981651332*Q1^2+8000000000000000000000000000000000*E11^3*B1^2+36*Q1*(23379999812960001496319988029439720749282234005742127954062976369*Q1^2-67957866123003737682636765205572000000000000000000000000000000000*E11^3*B1^2)^(1/2))*B1)^(1/3)-20000000000000000000000/3*E11^2*B1/((-5504587155963302752293577981651332*Q1^2+8000000000000000000000000000000000*E11^3*B1^2+36*Q1*(23379999812960001496319988029439720749282234005742127954062976369*Q1^2-67957866123003737682636765205572000000000000000000000000000000000*E11^3*B1^2)^(1/2))*B1)^(1/3));
                y10=[-100000,-100000,-100000];
                %plot(dz(l),y2c,'ro');
                if imag(y101)/abs(y101) < 0.000001 && real(y101) > 0
                    y10(1)=real(y101);
                end
                if imag(y102)/abs(y102) < 0.000001 && real(y102) > 0
                    y10(2)=real(y102);
                end
                if imag(y103)/abs(y103) < 0.000001 && real(y103) > 0
                    y10(3)=real(y103);
                end
                y11=fliplr(sort(y10));                    
                if y11(1)>0
                    plot(dz(l),y11(1),'m+');
                    y1sub(l)=y11(1);
                end
                if y11(2) >0
                    plot(dz(l),y11(2),'mx');
                    y1sup(l)=y11(2);                    
                end
                plot(dz(l),E11,'mo');
                E11up(l)=E11;
                E2=E11-dz2;
                y201=1/600000000000/B2*((-5504587155963302752293577981651332*Q2^2+8000000000000000000000000000000000*E2^3*B2^2+36*Q2*(23379999812960001496319988029439720749282234005742127954062976369*Q2^2-67957866123003737682636765205572000000000000000000000000000000000*E2^3*B2^2)^(1/2))*B2)^(1/3)+200000000000/3*E2^2*B2/((-5504587155963302752293577981651332*Q2^2+8000000000000000000000000000000000*E2^3*B2^2+36*Q2*(23379999812960001496319988029439720749282234005742127954062976369*Q2^2-67957866123003737682636765205572000000000000000000000000000000000*E2^3*B2^2)^(1/2))*B2)^(1/3)+1/3*E2;
                y202=-1/1200000000000/B2*((-5504587155963302752293577981651332*Q2^2+8000000000000000000000000000000000*E2^3*B2^2+36*Q2*(23379999812960001496319988029439720749282234005742127954062976369*Q2^2-67957866123003737682636765205572000000000000000000000000000000000*E2^3*B2^2)^(1/2))*B2)^(1/3)-100000000000/3*E2^2*B2/((-5504587155963302752293577981651332*Q2^2+8000000000000000000000000000000000*E2^3*B2^2+36*Q2*(23379999812960001496319988029439720749282234005742127954062976369*Q2^2-67957866123003737682636765205572000000000000000000000000000000000*E2^3*B2^2)^(1/2))*B2)^(1/3)+1/3*E2+1/200000000000*i*3^(1/2)*(1/6/B2*((-5504587155963302752293577981651332*Q2^2+8000000000000000000000000000000000*E2^3*B2^2+36*Q2*(23379999812960001496319988029439720749282234005742127954062976369*Q2^2-67957866123003737682636765205572000000000000000000000000000000000*E2^3*B2^2)^(1/2))*B2)^(1/3)-20000000000000000000000/3*E2^2*B2/((-5504587155963302752293577981651332*Q2^2+8000000000000000000000000000000000*E2^3*B2^2+36*Q2*(23379999812960001496319988029439720749282234005742127954062976369*Q2^2-67957866123003737682636765205572000000000000000000000000000000000*E2^3*B2^2)^(1/2))*B2)^(1/3));
                y203=-1/1200000000000/B2*((-5504587155963302752293577981651332*Q2^2+8000000000000000000000000000000000*E2^3*B2^2+36*Q2*(23379999812960001496319988029439720749282234005742127954062976369*Q2^2-67957866123003737682636765205572000000000000000000000000000000000*E2^3*B2^2)^(1/2))*B2)^(1/3)-100000000000/3*E2^2*B2/((-5504587155963302752293577981651332*Q2^2+8000000000000000000000000000000000*E2^3*B2^2+36*Q2*(23379999812960001496319988029439720749282234005742127954062976369*Q2^2-67957866123003737682636765205572000000000000000000000000000000000*E2^3*B2^2)^(1/2))*B2)^(1/3)+1/3*E2-1/200000000000*i*3^(1/2)*(1/6/B2*((-5504587155963302752293577981651332*Q2^2+8000000000000000000000000000000000*E2^3*B2^2+36*Q2*(23379999812960001496319988029439720749282234005742127954062976369*Q2^2-67957866123003737682636765205572000000000000000000000000000000000*E2^3*B2^2)^(1/2))*B2)^(1/3)-20000000000000000000000/3*E2^2*B2/((-5504587155963302752293577981651332*Q2^2+8000000000000000000000000000000000*E2^3*B2^2+36*Q2*(23379999812960001496319988029439720749282234005742127954062976369*Q2^2-67957866123003737682636765205572000000000000000000000000000000000*E2^3*B2^2)^(1/2))*B2)^(1/3));

                y20=[-100000,-100000,-100000];
                if imag(y201)/abs(y201) < 0.000001 && real(y201) > 0
                    y20(1)=real(y201);
                end
                if imag(y202)/abs(y202) < 0.000001 && real(y202) > 0
                    y20(2)=real(y202);
                end
                if imag(y203)/abs(y203) < 0.000001 && real(y203) > 0
                    y20(3)=real(y203);
                end                
                y2=fliplr(sort(y20));

                if dz2 > dzm
                   disp(dzm);
                end
                plot(dz(l),E2,'bo');
                E2down(l)=E2;
                if y2(1)>0
                    plot(dz(l),y2(1),'b+');
                    y2sub(l)=y2(1);
                end
                if y2(2) >0
                    plot(dz(l),y2(2),'bx');
                    y2sup(l)=y2(2);
                end                

                %plot(dz,(dz>-10000000)*yc,'b-'); 
                %plot(dz,(dz>-10000000)*E1,'y-');
                %plot(dz,(dz>-10000000)*E2c,'m-');
                %plot(dzm,0,'cd');
                %plot(dz(l),y2c,'m.');*/
                hold off
                figure(2), plot(dz,E11up,'m');
                hold on;
                figure(2), plot(dz,E2down,'c');
                figure(2), plot(dz,y1sub,'g');
                figure(2), plot(dz,y2sub,'r');
                figure(2), plot(dz,y1sup,'b');                
                figure(2), plot(dz,y2sup,'k');
            end
            pause(1);
        end
    end
end

                                    
                                

