function [sc_r,Cs_r,sc_w,Cs_w]=scoord0(theta_s,theta_b,Tcline,N)
% SCOORD0 returns just the functional form of the s-coordinate - no z values
% usage: [sc_r,Cs_r,sc_w,Cs_w]=scoord0(theta_s,theta_b,Tcline,N)

if nargin==1,
    theta_b = theta_s(2);
    Tcline  =theta_s(3);
    N       = theta_s(4);
    theta_s = theta_s(1);
end
    
    % code lifted from hernan's scoord3.m
    c1=1.0;
    c2=2.0;
    p5=0.5;
    Np=N+1;
    ds=1.0/N;
 
% rho points
    Nlev=N;
    lev=1:N;
    sc=-c1+(lev-p5).*ds;
    Ptheta=sinh(theta_s.*sc)./sinh(theta_s);
    Rtheta=tanh(theta_s.*(sc+p5))./(c2*tanh(p5*theta_s))-p5;
    Cs=(c1-theta_b).*Ptheta+theta_b.*Rtheta;
    sc_r = sc(:);
    Cs_r = Cs(:);    
    % w points
    Nlev=Np;
    lev=0:N;
    sc=-c1+lev.*ds;
    Ptheta=sinh(theta_s.*sc)./sinh(theta_s);
    Rtheta=tanh(theta_s.*(sc+p5))./(c2*tanh(p5*theta_s))-p5;
    Cs=(c1-theta_b).*Ptheta+theta_b.*Rtheta;
    sc_w = sc(:);
    Cs_w = Cs(:);