clear all;

n=[0.035];
s0=.0005;

%x=[0     -1772.60033     -3864.44619     -6436.35608     -8380.30788     -11867.4542];
%x=[0:-100:-10000];
x=[0:-200:-10000];

%y=[1.69 1.8 2 2.2 2.4 2.65 2.8 2.9 2.94 2.96];
%y=[4.5 3.9 3.4 3.1 3.03];
%yinv2=1./(y.^2);

q=3.8343;
m=0;
hold off;

    x1=0;
    y1=4.5;
    y3=y1;
    plot(x1,y1,'ro');
    g=9.81;
    
for x2=x(2:end)
    delx=x2-x1;
    v1=q/y1;

    sfk1=q^2*n^2/y1^(10/3);
    dydx1= (s0-sfk1)/(1-q^2/g/y1^3);
    k1=delx*dydx1;
    
    sfk2=q^2*n^2/(y1+k1/2)^(10/3);
    dydx2= (s0-sfk2)/(1-q^2/g/(y1+k1/2)^3);
    k2=delx*dydx2;
    
    sfk3=q^2*n^2/(y1+k2/2)^(10/3);
    dydx3= (s0-sfk3)/(1-q^2/g/(y1+k2/2)^3);
    k3=delx*dydx3;

    sfk4=q^2*n^2/(y1+k3)^(10/3);
    dydx4= (s0-sfk4)/(1-q^2/g/(y1+k3)^3);
    k4=delx*dydx4;
    
    y2=y1+1/6*(k1+2*k2+2*k3+k4);
   y3=[y3 y2];
    %str1=sprintf('x1=%d      x2=%d     delx=%d      so=%d             sfk1=%d     k1=%d             sfk2=%d      k2=%d            sfk3=%d      k3=%d         sfk3=%d    k4=%d        y1=%d      y2 = %d' ,x1,x2,delx,s0,sfk1,k1,sfk2,k2,sfk3,k3,sfk4, k4,y1,y2);
    str1=sprintf('x=%d       y = %d' ,x2,y2);
    disp(str1);
    hold on;
    %plot(x2,y2,'ro');
     y1=y2;
     x1=x2;
    pause(.1)
end
    plot(x,y3,'r');

    
    
y=[4.5:-.005:3.0]
q=3.8343;
m=0;

R=y;

sf=q^2*n^2./R.^(10/3)

sfa=[sf(1:end-1)+sf(2:end)]/2;

s0msfa=s0-sfa;
v=q./y;
E=y+v.^2/(2*9.81);
dE=E(2:end)-E(1:end-1);

x0=0;

delx=[0 dE./s0msfa];
x=cumsum(delx);
str1=sprintf('y = %s             \nx=%s\n',num2str(y),num2str(x));
disp(str1);
%hold off;
%plot(x,y,'bo');
hold on;
plot(x,y,'b');
pause(1)


x=[0:200:10000];
    x1=0;
    y1=4.5;
    y3=y1;
    plot(x1,y1,'yo');
    g=9.81;
   str1=sprintf('x = %d      y=%d',x1,y1) ;
   disp(str1);
for x2=x(2:end)
    delx=x2-x1;
    sf1=q^2*n^2/y1^(10/3);
    v1=q/y1;
    rhs=    y1+v1^2/(2*g)-s0*delx+0.5*sf1*delx;
    EqnString=sprintf('y2+%d^2/(2*9.81*y2^2)-0.5*(%d^2*%d^2/y2^(10/3))*(%d)=%d',q,n,q,delx,rhs);
    
    %instr1=sprintf('Q2=%d,Q1=%d,B=%d,g=9.81,m=%d,n=%d,A2=(B+m*y2)*y2,P2=(B+2*sqrt(m^2+1)*y2),R2=A2/P2,A1=(B+m*y1)*y1,P1=(B+2*sqrt(m^2+1)*y1),R1=A1/P1,s0=%d,delx=%d,y1=%d',Q,Q,B,m,n,s0,delx,y1)
    %ANS1=solve('y2+Q2^2/(2*g*A2^2)+0.5*(n^2*Q2^2/(A2^2*R2^(4/3)))*delx=s0*delx+y1+Q1^2/(2*g*A1^2)-0.5*(n^2*Q1^2/(A1^2*R1^(4/3)))*delx',instr1);
    ANS1=solve(EqnString);
   y2=ANS1(find(~(imag(ANS1)~=0)>0));
   str1=sprintf('x = %d      y=%d',x2,double(y2));
      disp(str1);
   x1=x2;
    y1=double(y2);
    y3=[y3 y1];
    hold on;
    %plot(-x2,double(y2),'go');
    pause(.1)
end
    plot(-x,y3,'g');
