clear all;

%Radius of the particle range
RawRp = [1,25,50,70,100,500]*1e-9;

%Contact angle range
thetapRaw = [50, 25, 50,  50];   
thetasRaw = [10, 5,  2.5, 5 ];


for j=1:size(thetapRaw,2)
   
    for i=1:size(RawRp,2)

        Rp = RawRp(i);

        andri = Rp/100;
        thetap=thetapRaw(j)*pi/180;
        thetas=thetasRaw(j)*pi/180;

        z00=1e-9;
        for p=1:101
            z=Rp-p*andri+andri; 
            Z(p)=z;
            Z0(p)=fsolve(@(z0) elevation(z0,Rp,thetap,thetas),z00);
            z0=Z0(p);

            % r(z)
            a=-asinh(cot(thetas));
            beta=(Rp-z0)/sqrt(Rp^2-(Rp-z0)^2);
            b=1/z0*(1/a*asinh((beta*cos(thetap)-sin(thetap))/(beta*sin(thetap)+cos(thetap)))-1);
            r=1/(a*b)*cosh(a*(b*z+1));
            R(p)=r;

        end
        maxR(i) = max(R)+Rp;

        if (Rp>20e-9)&&(Rp<200e-9)
            if (thetap==50*pi/180)&&((thetas==2.5*pi/180) ||(thetas==10*pi/180))
                figure(101); hold on;
                plot(R*1e+6,Z*1e6);
                xlabel('\itr\rm [\mum]');
                ylabel('\itz\rm [\mum]');
            end
        end

    end
    figure(102); hold on;
    plot(RawRp*1e+6,maxR*1e6);
    xlabel('\itRp \rm [\mum]');
    ylabel('\itz\rm [\mum]');
    drawnow
end