% Script splinmat made from PROGRAM SPLINRUN in Pascal
% Described in Malacic V.: Ecological Modeling, 138, 173-191, 2001 and 
% Malacic V., 1998. Geofizikalno-ekološki pristop k disperziji odplak piranskega izpusta. 
% Fazno porocilo 2, Nacionalni inštitut za biologijo, Morska biološka postaja Piran, 43 str.
      NN = 500; % maximum dimension of arrays, also number of steps in Runge-Kutta}
      kmax = 400;%200}; % max. number of results, kount < kmax}
      NV = 4; % Num. of  variables }

% An Introduction to Water Quality Modelling, A.James (Ed.), Wiley, 1984, 234 pp., 
% page 155: Lambda = 1.16;  Fischer et al., 1979, Mixing in Inland Waters and Coastal Waters,
% Academic Press, page 370: Lambda = 1.2
      Lambda = 1.14; % Lambda = ratio between the radius at which the density 
      % difference (ambient density minus core density of the buoyant jet) falls 
      % to 1/e of the core density difference, and the radius at which the 
      % velocity falls to 1/e of the core velocity in a buoyant slice. In both, the density 
      % difference profile and the velocity profile across the buoyant jet the Gaussian profile is supposed. 
      % for Lamba =1.14 the best numerical results were obtained when compared with semiempirical solutions.
      ystart = zeros(NV,1);% starting values of variables; glnarray;
      % yp = size NN x NV = results, four output variables, see below 'initial data'
      time = zeros(NN,1);
% begin % MAIN program  }
%**************************************************}
   varyAlpha = input('Do you want to vary Alpha (the enrtainment coefficient) as a function of Ri number? Y/N [N]:','s');
   % Fischer et al., 1979, Mixing in Inland Waters and Coastal Waters, Academic Press
   if isempty(varyAlpha)
      varyAlpha = 'N';
   end
   varyAlpha = upper(varyAlpha);
   
   DepthDifstr = '10.1'; % DepthDif = 10.1; %36.0;}%20.7;} % Depth of diffusser
   Diameterstr = '0.1';%'0.6';%0.1;}%3.0;;} % [m] diameter of the orifice } 0.6 m = no difussor
   u0str = '12.7';%'3.54';%12.7;}%20.0} %1.27;} %0.5;} % [m/s] 3.54 m/s = no difussor, 1000 l/s
   theta0str =  '0';%PI*(90.0)/180.} PI*(0.0)/180.;   
   prompt = {'Enter diffusor depth (m):','Enter orifice diameter (m):',...
             'Enter initial speed from the orifice (m/s):',...
             'Enter initial jet inclination above the horizontal plane (degrees):'};
   dlgtitle = 'Input params';  
   dims = [1 30; 1 35; 1 45; 1 70];% [1 70]
   definput = {DepthDifstr,Diameterstr,u0str,theta0str};% dif. depth(m), orifice diam.(m), u0(m/s), theta(deg)
   replies = inputdlg(prompt, dlgtitle, dims,definput);
   initstrngs = string(replies);
   DepthDif = str2double(initstrngs(1));
   Diameter = str2double(initstrngs(2));
   u0 = str2double(initstrngs(3));
   theta0 = str2double(initstrngs(4));
   theta0 = theta0*pi/180;
   SS0 = 6.2*Diameter;
   s1 = 0.001; % guessed first stepsize }
   smin = 1.0E-6; % minimum step size = 1 micron, is 'small enough'}
   SSmax = 100.0;% not known final distance along the jet }
   dxsav = 0.1;%recorded points are separated by ds > 0.1 }
   b0 = Diameter/2^0.5; % Water Quality Model, p.155 }
   bmax = 100,0; % maximum radius of the buoyant jet - chosen arbitrary to close the 
   % program if radius of the jet becomes too large before the 'final'
   % depth (or sea surface) is achieved
   [NumDat,HeightArr,GamArr] = ReadData(DepthDif);

   %***** GET DENSITY AT INITIAL HEIGHT gam0 BY INTERPOLATION FROM THE DENSITY PROFILE *****}
   Height0 = eps; % starting height above the diffusser }
   tempklo0 = find((HeightArr(1:NumDat) <= Height0),1,'last')
   % get coefficients of polynoms in each interval withing points HeightArr
  
   YesMakima = input('Do you want to apply the Makima spline (Matlab R2019b) Y/N. Default [N]:','s');
   if isempty(YesMakima)
      YesMakima = 'N';
   end
   YesMakima = upper(YesMakima);

   % get the interpolated value of density gam0 and its vertical derivative 
   % dgamdz at the diffuser height Height0 
   if YesMakima == 'Y'
      ppmakima = makima(HeightArr,GamArr);
      %https://blogs.mathworks.com/cleve/2019/04/29/makima-piecewise-cubic-interpolation/
      [gam0, dgamdz0]=ppvalder(ppmakima,Height0);
   else
      ppsplin = spline(HeightArr,[0 GamArr' 0]);
      % [0 GamArr 0] added zero first derivatives in the first and the last data point     
      [gam0, dgamdz0]=ppvalder(ppsplin,Height0);
   end   
   
   % 
   fprintf(1,'HeightArr(%d)= %6.2f Height0 = %6.2f gam0= %7.3f dgamdz0= %7.3f\r\n',...
           tempklo0,HeightArr(tempklo0),Height0, gam0, dgamdz0);     
   %************   INITIAL DATA  *******************}
   % SS0 starting distance [m] along the jet }
   % ystart[1] = velocity along the jet axis }
   % ystart[2] = half-width of a buoyant jet }
   % ystart[3] = angle of jet axis with respect to XX }
   % ystart[4] = delta_rho = rho_ambient - rho_jet }
   ystart(1) = u0;
   ystart(2) = b0;
   ystart(3) = theta0;
   ystart(4) = gam0*(1.0 + (Lambda)^2)/(2.0*(Lambda)^2);
   
   % suppose a sewage with the density of 1000 kg/m^3,  }
   % delta_Rho0 = gam0, also 'Water Quality Model.',1984, p.154 }
   Rho0 = 1000.0 + gam0;%Rho0 = 995.8115 + gam0;} %T=30 st.C, S=0, p=36 dbar}
   

   %*********   NAMES OF OUTPUT FILES ************}
   prompt = {'Enter the output file name for buoyant jet drawing:'};
   dlgtitle = 'output data for drawing';
   dims = [1 35];
   defoutdrw = {'outputdraw.txt'};
   answer = inputdlg(prompt,dlgtitle,dims,defoutdrw);   
   fiddrw = fopen(string(answer),'w+');
   
   prompt = {'Enter the output file name for key output data:'};
   dlgtitle = 'output key data';
   dims = [1 35];
   defoutdrw = {'outputkey.txt'};
   answer = inputdlg(prompt,dlgtitle,dims,defoutdrw);   
   fidkey = fopen(string(answer),'w+');
   
   prompt = {'Enter the output file name for comparing num and analytical results:'};
   dlgtitle = 'output data for comparison';
   dims = [1 35];
   defoutdrw = {'outputcom.txt'};
   answer = inputdlg(prompt,dlgtitle,dims,defoutdrw);   
   fidcom = fopen(string(answer),'w+');
  
%(***   CALL THE BEST FAECALI ENGINE:   *****)
%**************************************************}
   [nok,nbad,kount,xp,yp,xxp,zzp,RiArr,AlphArr] = FAECALI(ystart,...
                  NV,kmax,...
                  HeightArr,GamArr,...
                  varyAlpha,DepthDif,...
                  SS0,SSmax,s1,smin,bmax,dxsav,...
                  Lambda,Diameter,Height0,theta0,YesMakima);
   
   fprintf(1,'success. steps (nok)= %d  bad steps (nbad)= %d stored interm. values (kount) = %d\r\n',...
           nok,nbad,kount);
   xp = xp(1:kount); xxp = xxp(1:kount); zzp = zzp(1:kount); % shorten arrays to results
   RiArr = RiArr(1:kount); AlphArr = AlphArr(1:kount);
   yp = yp(1:kount,:);
   for i =2 : kount % kount 
      time(i) = time(i-1) + ...
      (xp(i)-xp(i-1))*2.0/(yp((i),1)+yp(i-1,1));
   end
   % time = time(1:kount); 
   time = time(1:kount); 
   
%(***   Call the STRATIMIDSURF function to write down results   *****)
%**************************************************} 
   [Sm,SSstra,ZZstra,minkount] = STRATIMIDSURF(xp,yp,time,xxp,zzp,NV,kount,...
                        Rho0,Height0,DepthDif,RiArr,AlphArr,Lambda,Diameter,u0,...
                        HeightArr,GamArr,fiddrw,fidkey,fidcom,YesMakima);
   % Sm = scalar, the last value of numerical dilution. SSstra = calculated dilution
   % from analytic expression for a linear vertical density profile. ZZstra = max height of sewage, 
   % calculated analytiucally from a linear vertical density profile.
   % if zzp(kount) < zmax = DepthDif 
   %      then zzp(del_rho=0) IS WITHIN the WATER COLUMN
   % else at the surface 
   fprintf(1,'After STRATIMIDSURF kount=%d\r\n', kount);
   fclose(fiddrw);
   fclose(fidkey);
   fclose(fidcom);
 %*************************** put data on screen and plot the data  
   disp('  time(s)    xp(m)    yp1(m/s)    yp2(m)   yp3(rad)  yp4(kg/m3)  xxp(m)   zzp(m)');
   [time(1:10) xp(1:10) yp(1:10,:) xxp(1:10) zzp(1:10)]
   disp('  time(s)    xp(m)    yp1(m/s)    yp2(m)   yp3(rad)  yp4(kg/m3)  xxp(m)   zzp(m)');
   [time(kount-1:kount) xp(kount-1:kount) yp(kount-1:kount,:) xxp(kount-1:kount) zzp(kount-1:kount)]
   figure(1)
   plot(time(1:kount),yp(1:kount,1),'k',time(1:kount),yp(1:kount,2),'b',...
        time(1:kount),yp(1:kount,3).*(180/pi),'g',time(1:kount),yp(1:kount,4),'r','LineWidth',2 )
   xlabel('{\it{time}}(s)');
   legend('\it{u}(m/s)','\it{b}(m)','\theta (°)','{\Delta\rho}(kg/m^3)');
   grid on;
   figure(2)
   xxpbr = xxp+yp(:,2).*sin(yp(:,3)); zzpbr = zzp-yp(:,2).*cos(yp(:,3)); % br = bottom right
   indbr = find(zzpbr <= DepthDif,1,'last'); xxpbr = xxpbr(1:indbr); zzpbr = zzpbr(1:indbr); 
   xxptl = xxp-yp(:,2).*sin(yp(:,3)); zzptl =  zzp+yp(:,2).*cos(yp(:,3)); % tl = top left
   indtl= find(zzptl <= DepthDif,1,'last'); xxptl = xxptl(1:indtl); zzptl = zzptl(1:indtl); 
   indc = find(zzp <= DepthDif,1,'last'); % c = central
   xxpc = xxp(1:indc); zzpc = zzp(1:indc);
   plot(xxpc,zzpc,'r',xxpbr,zzpbr,'b:',xxptl,zzptl,'b:','LineWidth',2);
   xlabel('{\it{x}}(m)');ylabel('{\it{z}}(m)');
   grid on; axis equal;
   figure(3)
   plot(time,zzp,'LineWidth',2)
   xlabel('{\it{time}}(s)');ylabel('{\it{z}}(m)');
   grid on;
   
   
function [NumDat,HeightArr,GamArr] = ReadData(DepthDif)
%--------------------------------------------------}
   filelist = dir;
   filelist = filelist(endsWith({filelist.name}, 'txt'));
   inpfn = {filelist.name};
   [indx] = listdlg('PromptString',{'Select the input file.',...
       'Only one file with columns of depth and density',''},...
       'SelectionMode','single','ListString',inpfn);
   finp = fopen(filelist(indx).name,'r');
   % formatSpec = '%f %f %f %f';
   formatSpec = '%f %f %*f %*f'; % skip last two real values (temperature and salinity)
   % sizeInp = [4 Inf];
   sizeInp = [2 Inf];
   InpDat = fscanf(finp,formatSpec,sizeInp)';
   NumDat = length(InpDat(:,1)); % get the number of data pairs
   
   HeightArr = DepthDif - flip(InpDat(:,1)); % height above the diffusser
   GamArr = flip(InpDat(:,2));
   fclose(finp);
end % ReadData }


function [Ri, Alpha] = MakeRiAlpha(y, Alphj, Alphp, Lambda, Rip, Rho0, GG, varyAlpha)
%**************************************************}
  SqrRi = 4.0*(Lambda)^2*(2.0*pi)^0.5*GG*y(2)*y(4)/...
          (Rho0*(1.0+(Lambda)^2)*(y(1))^2);
  if (y(2)*y(4) > 0) % then
    Ri = (SqrRi)^0.5;
  else
    Ri = -(-SqrRi)^0.5;
  end  
  % if UpCase(% varyAlpha) = 'Y' % then
  if varyAlpha == 'Y' % added: = 'Y'
    if SqrRi < -200.0 % then
      Alpha= 0.0;
    else
      Alpha= Alphj*exp(log(Alphp/Alphj)*SqrRi/(Rip)^2); 
      % Fischer et al., 1979, Mixing in Inland Waters and Coastal Waters, 
      % Academic Press, page 371, expression (9.144) in a Chapter 9,
      % 'Turbulent jets and plumes'
    end  
  else
    Alpha = Alphp; % Alpha does not vary, does not depend on Ri (Richardson) number
  end  
end % MakeRiAlpha }


function dydx = DERIVS(~,y,GG,Rho0,dgamdz,Alpha,Lambda) % ~ stands for 'x', or 's' - 
%*** the independent argument, zthe length of trajectory, which is not present in dydx }
   % y[1] = velocity along the jet axis }
   % y[2] = half-width of a buoyant jet }
   % y[3] = angle of jet axis with respect to XX }
   % y[4] = delta_rho = rho_ambient - rho_jet }
   dydx(1) = 2.0*GG*(Lambda)^2*y(4)*sin(y(3))...
             /(y(1)*Rho0) - 2.0*y(1)*Alpha/y(2);
   dydx(2) = 2.0*Alpha-(y(2)*GG*(Lambda)^2*y(4)*sin(y(3))...
             /((y(1)^2)*Rho0));
   dydx(3) = 2.0*GG*(Lambda)^2*y(4)*cos(y(3))...
             /((y(1)^2)*Rho0);
   dydx(4) = (1.0+(Lambda)^2)*sin(y(3))*dgamdz...
             /((Lambda)^2) - 2.0*Alpha*y(4)/y(2);
end

function yout = rk4(y,dydx,NV,x,h,GG,Rho0,dgamdz,Alpha,Lambda)
%**************************************************}
   yt = zeros(NV,1)'; yout = yt;
   hh = h*0.5;
   h6 = h/6.0;
   xh = x+hh;
   for i = 1 :NV % do % begin
      yt(i) = y(i)+hh*dydx(i);
   end
   dyt = DERIVS(xh,yt,GG,Rho0,dgamdz,Alpha,Lambda);
   for i = 1 : NV % do % begin
      yt(i) = y(i)+hh*dyt(i);
   end
   dym = DERIVS(xh,yt,GG,Rho0,dgamdz,Alpha,Lambda);
   for i = 1 : NV % do % begin
      yt(i) = y(i)+h*dym(i);
      dym(i) = dyt(i)+dym(i);%  
   end
   dyt = DERIVS(x+h,yt,GG,Rho0,dgamdz,Alpha,Lambda);
   for i = 1 : NV % do % begin
      yout(i) = y(i)+h6*(dydx(i)+dyt(i)+2.0*dym(i));
   end
end

function [x,y, hdid, hnext] = RKQC(ystart,dydx,NV,xsav,htry,EPS,yscal,...
                                 GG,Rho0,dgamdz,Alpha,Lambda)
%**************************************************}
% provide for the call the derivatives dydx at location x, given
% both x and the function values y. They must also define
% added ystart in the call on 8 Nov. 2019
   pgrow=-0.20;
   pshrnk=-0.25;
   fcor=0.06666666666666666666;   %(* 1.0/15.0 *)
   one=1.0;
   safety=0.9;
   errcon=6.0e-4;
   ysav = zeros(NV)'; dysav = ysav;
   for i = 1 : NV % do % begin
      % ysav(i) = y(i);
      ysav(i) = ystart(i);
      dysav(i) = dydx(i);
   end
   h = htry;
   errmax = 1.1; % added to start the while loop 
   x = xsav + 0.1; % added to start the while loop
   % while (errmax > one)
   while ((errmax > one) & (x-xsav)> EPS)
      hh = 0.5*h;
      % evaluate first Runge-Kutta for a step h/2:}
      ytemp = rk4(ysav,dysav,NV,xsav,hh,GG,Rho0,dgamdz,Alpha,Lambda);
      x = xsav+hh;
      dydx = DERIVS(x,ytemp,GG,Rho0,dgamdz,Alpha,Lambda);
      % second step starts at x=xsav+0.5*h. }
      % Evaluate second Runge-Kutta for a step h/2:}
      y = rk4(ytemp,dydx,NV,x,hh,GG,Rho0,dgamdz,Alpha,Lambda); % you have obtained y }
      x = xsav+h;
      if (abs(x-xsav) <= EPS) % then % begin
         input('stop in routine RKQC\r\n stepsize too small. Press Return to continue','s');
         % return;
      end
      % Evaluate Runge-Kutta for a full step h:}
      ytemp = rk4(ysav,dysav,NV,xsav,h,GG,Rho0,dgamdz,Alpha,Lambda);% ytemp by large step}
      errmax = 0.0;
      for i = 1 : NV % do % begin
         ytemp(i) = y(i)-ytemp(i);% difference in methods}
         temp = abs(ytemp(i)/yscal(i)); % relative error }
         if (errmax < temp) % then 
             errmax = temp;
         end    
      end
      errmax = errmax/EPS;
      % if (errmax > one) % then % begin
      h = safety*h*exp(pshrnk*log(errmax));
   end % while
   hdid = h;
   if (errmax > errcon) % then % begin
      hnext = safety*h*exp(pgrow*log(errmax));
   else % begin
      hnext = 4.0*h;
   end
   for i = 1 : NV
      y(i) = y(i)+ytemp(i)*fcor;
   end
end % RKQC }

function [nok,nbad,kount,xp,yp,xxp,zzp,RiArr,AlphArr] = FAECALI(ystart,...
                  NV,kmax,...
                  HeightArr,GamArr,...
                  varyAlpha,DepthDif,...
                  x1,x2,s1,smin,bmax,dxsav,...
                  Lambda,Diameter,Height0,theta0,YesMakima)
%**************************************************}
   Zmax = DepthDif;
   GG = 9.806; % for the TS Gulf }
   EPS = smin; % accuracy }
   EPS_U = 1.0e-5;
   Alphj = 0.0535; % Fischer et al., 1979, page 371: }
   Alphp = 0.0833;
   Rip = 0.557;
   maxstp=10000;
   NN = 500;
   two=2.0;
   zero=0.0;
   tiny=1.0e-30;
   xp = zeros(NN,1); xxp = xp; zzp = xp; RiArr = xp; AlphArr = xp;
   yp = zeros(NN,NV);
   x = x1;
   xxp(1) = 6.2*Diameter*cos(theta0);% starting distance }
   zzp(1) = Height0+6.2*Diameter*sin(theta0);
   if (x2 > x1) % then 
      ds = abs(s1);
   else
       ds = -abs(s1);
   end
   nok = 0;
   nbad = 0;
   kount = 1;%1 0;
   y = ystart;
   yscal = zeros(1:NV);
   xp(1) = x1; % inserted on 15 September 98 } Nov 2019: x1 = SS0 = 3.72
   xsav = x-dxsav*two; % dxsv = 0.1
   yold = ystart;% 14 Nov 2019: to define the 'previous values' of solution 
   yp(1,:) = ystart;

   
   if YesMakima == 'Y'
      ppppspl = makima(HeightArr,GamArr');
      % https://blogs.mathworks.com/cleve/2019/04/29/makima-piecewise-cubic-interpolation/
   else
      ppppspl = spline(HeightArr,[0 GamArr' 0]); 
      % [0 GamArr' 0] added zero first derivatives in the first and the last data point
   end   
   
% (******     START BIG for LOOP:    ****************)
for nstp = 1 : maxstp % do % begin
   [gama, dgamdz]=ppvalder(ppppspl,zzp(kount));
   fprintf(1,'Start loop, nstp= %3d; kount= %3d; xp(%3d)= %6.4f; xxp(%3d)= %6.4f; zzp(%3d)= %6.4f; gama= %6.2f ; dgamdz= %10.4f\n',...
           nstp,kount,kount,xp(kount),kount,xxp(kount),kount,zzp(kount),gama,dgamdz);
   Rho = gama + 1000;
   [Ri, Alpha] = MakeRiAlpha(y,Alphj,Alphp,Lambda,Rip,Rho,GG,varyAlpha);
   fprintf(1,'x= %7.4f Ri= %6.2f ; Alpha= %6.2f\n',x,Ri,Alpha);
   fprintf(1,'y(1)= %10.5f ; y(2)= %10.5f ; y(3)= %10.5f ; y(4)= %10.5f\n',...
            y(1),y(2),y(3),y(4));
   dydx = DERIVS(x,y,GG,Rho,dgamdz,Alpha,Lambda);
   fprintf(1,'dydx(1)= %10.5f ; dydx(2)= %10.5f ; dydx(3)= %10.5f ; dydx(4)= %10.5f\n',...
           dydx(1),dydx(2),dydx(3),dydx(4));
   fprintf(1,'kount=%4d x=%7.4f xsav=%7.4f dxsav=%5.4f zzp(%4d)=%7.4f; Zmax=%6.2f\n\n',...
                kount,x,xsav,dxsav, kount,zzp(kount),Zmax);
   for i = 1 : NV% var % do % begin
  % Scaling for monitoring accuracy - modify if necessary }
         yscal(i) = abs(y(i))+abs(dydx(i)*ds)+tiny;
   end
   xsav = x;
   yold = y; 
   [x,y, dsdid, dsnext] = RKQC(yold,dydx,NV,xsav,ds,EPS,yscal,...
                               GG,Rho,dgamdz,Alpha,Lambda);% added ystart      
   ds = dsnext;                        
   %% nok = number of good steps, }
   %% nbad = number of bad (repaired) steps }
   %% kount = number of stored (written) values }
   if (dsdid - ds <= tiny) % then % begin
         nok = nok+1;
   else % begin
         nbad = nbad+1;
   end
   if ( (abs(x-xp(kount)) > abs(dxsav)) && (kount < kmax-1) && ...
         (((x-x2)*(x2-x1)) < zero)   && ...% x<x2=SSmax and x2=SSmax > x1=SS0
         ((y(2) <= bmax ) && (zzp(kount) > 0)) && ...
         ((y(3) >= 0) && (zzp(kount) > 0)) && ...
          (zzp(kount) < Zmax) )
         kount = kount+1;
         xp(kount) = x;
         for i = 1 : NV% var % do % begin
            yp(kount,i) = y(i); %
         end
         if kount > 1 % then % begin
            xxp(kount) = xxp(kount-1) + ...
               (xp(kount)-xp(kount-1))*cos(y(3));
            zzp(kount) = zzp(kount-1) + ...
               (xp(kount)-xp(kount-1))*sin(y(3));
         else % kount = 1
            xxp(kount) = (xp(kount))*cos(y(3));
            zzp(kount) = (xp(kount))*sin(y(3));
         end    
            RiArr(kount) = Ri;
            AlphArr(kount) = Alpha;
   end % if (abs(x-xsav) > abs(dxsav)) }
   % if stepsize can overshoot, decrease it }
   if (((x+ds-x2)*(x+ds-x1)) > zero) % then x+ds > x2 and x+ds>x1; x2 = SSmax = 100
      ds = x2-x;
   end     
   % Are we done?: if (((x-x2)*(x2-x1)) >= zero) THEN % begin
   % finish when y[1] = u < 0, NOT when theta (y[3]) < 0, 
   % as supposed in Water Quality Modelling (p. 162) }
   maxzzp = max(zzp(1:kount));
   if ( (((x-x2)*(x2-x1)) >= zero) || ...% x>x2=SSmax and x2=SSmax > x1=SS0 
      ((y(3) < 0.0) && (zzp(kount) > 0)) || ...
      ((y(1) < EPS_U) && (zzp(kount) > 0)) || ...
      ((y(2) > bmax ) && (zzp(kount) > 0)) || ...
      ((zzp(kount) < maxzzp) || (zzp(kount) >= Zmax)) ) 
         fprintf(1,'Closing LONG IF zzp(%4d)=%7.4f; Zmax=%6.2f\n',kount,zzp(kount),Zmax);
         fprintf(1,'Before the RETURN (=EXIT from Faecali) command;  kount= %6d ',kount); 
         prompt = 'Press Return\n\n'; str = input(prompt,'s');
         kount = kount+1;
         xp(kount) = x;
         for i = 1 : NV% var % do % begin
            yp(kount,i) = y(i); % 
         end
         if kount > 1 % then % begin
            xxp(kount) = xxp(kount-1) + ...
               (xp(kount)-xp(kount-1))*cos(y(3));
            zzp(kount) = zzp(kount-1) + ...
               (xp(kount)-xp(kount-1))*sin(y(3));
         else % kount = 1
            xxp(kount) = (xp(kount))*cos(y(3));
            zzp(kount) = (xp(kount))*sin(y(3));
         end    
            RiArr(kount) = Ri;
            AlphArr(kount) = Alpha;         
         return; % !!!! in Fortran (and Pascal) GOTO 99 % normal exit from function}
   end % if (()) >= 0 OR ... OR ... }
   if (abs(dsnext) < smin) % then % begin
         prompt = 'Pause in routine FAECALI, stepsize too small. Press Return!';
         str = input(prompt,'s');
         return; % exit when the step is too small
   end
end % for nstp }
% (***********   end BIG for LOOP *******************)
end % FAECALI } 


function [Sm,SSstra,ZZstra,minkount] = STRATIMIDSURF(xp,yp,time,xxp,zzp,NV,kount,...
                              Rho0,Height0,DepthDif,RiArr,AlphArr,Lambda,Diameter,U0,...
                              HeightArr,GamArr,fiddrw,fidkey,fidcom, YesMakima)
%*********************************************************************************************}
   % prompt = 'INSIDE STRATIMIDSURF. Press Return!\n'; str = input(prompt,'s');
   disp('START STRATIMIDSURF');
   GG = 9.81;
   fprintf(fidkey,'%s\r\n','  zzp   Sm    y(1)    y(2)   y(3)  y(4)');  
   for jj = 1 : NV % do
      fprintf(1,'jj=%3d; yp(%3d,%3d)= %6.3f\n',jj,kount,jj,yp(kount,jj));
   end 
   del_dens=diff(yp(:,4)); % look for differences of densities to find the first with del_dens = 0
   % because you need different densities to apply the spline function - which later you did not apply!
   minkount = find(del_dens >= 0.0,1,'first');
   fprintf(1,'In STRATIMIDDLE; minkount=%4d\n',minkount);  
   if zzp(kount) < DepthDif % then % begin 
      prompt = 'REMAINED BELOW. Oki-Doki. Press Return!\n'; str = input(prompt,'s');
   else % tempklo >= kount % begin       
      disp('GONE TO THE SURFACE, Bye-bye.');
      fprintf(1,'kount=%4d; xp(kount)=%8.4f',kount,xp(kount));
 %%%   time  xp    yp1     yp2   yp3     yp4    xxp    zzp   xxp+r   zzp-r xxp-r zzp+r Sm    Ri    Alfa    
      size(yp)
      for jj = 1 : NV % do
          fprintf(1,'jj=%3d; yp(%3d,%3d)= %6.3f\n',jj,kount,jj,yp(kount,jj));
      end   
      fprintf(1,'kount=%4d; xxp(kount)= %5.2f; zzp(kount)= %5.2f; DepthDif= %5.2f\n',...
                  kount,xxp(kount), zzp(kount), DepthDif);
      prompt = 'before end if tempklo < kount; tempklo <= kount! Press Return!\n'; str = input(prompt,'s');
     
   end % if zzp(kount) < DepthDif  else ...
   %****************************************************************************************************}
   fprintf(fiddrw,'%s\r\n','time(s) xp(m)  yp1     yp2     yp3    yp4     xxp   zzp  xxpb  zzpb  xxpt  zzpt   Sm     Ri   Alpha');
   for i = 1 : kount % do % begin
      fprintf(fiddrw,' %5.2f %5.2f',time(i),xp(i));
      for jj = 1 : NV % do      % begin
          fprintf(fiddrw,' %6.3f ',yp(i,jj));
      end
      fprintf(1,'xxp(%3d)%5.2f zzp(%3d)=%5.2f\n',i,xxp(i),i,zzp(i));
      % Sm - numerical dilution along the trajectory }
      Sm = yp(i,1)*(yp(i,2)^2)*4.0*(Lambda)^2/...
           (U0*(Diameter)^2*(1.0+(Lambda)^2));
      fprintf(fiddrw,'%5.2f %5.2f %5.2f %5.2f %5.2f %5.2f %5.2f %6.3f %6.4f\r\n',...
              xxp(i),zzp(i),xxp(i)+yp(i,2)*sin(yp(i,3)),zzp(i)-yp(i,2)*cos(yp(i,3)),...
              xxp(i)-yp(i,2)*sin(yp(i,3)), zzp(i)+yp(i,2)*cos(yp(i,3)),...
              Sm,RiArr(i),AlphArr(i));
   end %i = 1:kount
   
   Sm = yp(1,1)*(yp(1,2))^2*4.0*(Lambda)^2/...
        (U0*(Diameter)^2*(1.0+(Lambda)^2)); % Sm = initial dilution here
   fprintf(fidkey,'%5.2f %5.2f %6.4f %6.4f  %6.4f %6.4f\r\n',...
           zzp(1),Sm,yp(1,1),yp(1,2),yp(1,3),yp(1,4));
   fprintf(fidcom,'1   zzp( 1 )[m]=%7.2f; Sm=%7.2f; yp( 1 ,1)[m/s]=%7.2f yp( 1 ,4)[m]=%7.2f\r\n',...
           zzp(1),Sm,yp(1,1),yp(1,4));   
   fprintf(1,'pred Sm, kount= %4d; zzp=%5.2f yp1=%5.2f; yp2=%6.4f; yp3=%6.4f; yp4=%6.4f\r\n',...
           kount,zzp(kount),yp(kount,1),yp(kount,2),yp(kount,3),yp(kount,4));       
   Sm = yp(kount,1)*(yp(kount,2))^2*4.0*(Lambda)^2/...
        (U0*(Diameter)^2*(1.0+(Lambda)^2)); % Sm = final dilution here
   fprintf(fidkey,'%5.2f %5.2f  %6.4f %6.4f %6.4f %6.4f\r\n',...
           zzp(kount),Sm,yp(kount,1),yp(kount,2),yp(kount,3),yp(kount,4));
   fprintf(fidcom,'%3d zzp(%3d)[m]=%7.2f; Sm=%7.2f; yp(%3d,1)[m/s]=%7.2f yp(%3d,4)[m]=%7.2f\r\n',...
           kount,kount,zzp(kount),Sm,kount,yp(kount,2),kount,yp(kount,4));
   
   %** <dgamdz> from diffusser to max plume height }
   gam0 = Rho0 - 1000;
   if YesMakima == 'Y'
      ppmakima = makima(HeightArr,GamArr);
      gama = ppval(ppmakima,DepthDif); 
   else
      ppsplin = spline(HeightArr,[0 GamArr' 0]);
      % [0 GamArr 0] added zero first derivatives in the first and the last data point     
      gama = ppval(ppsplin,DepthDif);
   end   

   % gama =ppval(ppmakima,DepthDif);
   dgamdz = (gama-gam0)/(DepthDif-Height0);                   
   fprintf(1,'<dgamdz>= %6.2f\r\n', dgamdz);
   %** call analytical estimation of ZZstra **** and SSstra - fixed dgamdz **}
   [ZZstra,SSstra] = ANALYTIC(Diameter,dgamdz,Rho0,U0,GG,DepthDif);
   fprintf(1,'ZZStra=%7.3f; maxzzp=%7.3f; SSstra= %5.2f; Sm= %5.2\r\n',... 
              ZZstra,zzp(kount),SSstra,Sm);
   fprintf(fidcom,'ZZStra=%7.3f; maxzzp=%7.3f; SSstra= %5.2f; Sm= %5.2f\r\n',... 
              ZZstra,zzp(kount),SSstra,Sm);
end % STRATIMIDSURF }


function [ZZstra,SSstra] = ANALYTIC(diameter,dgamdz,Rho0,u0,GG,DepthDif)
%************ calculate dilution supposing the linear density profile 
%************ from the height of diffusser to the max achieved height
% const
      konst1 = 3.98;
      konst2 = 0.071;
      konst3 = 0.089;
   delRho0 = Rho0 - 1000.0;
   QQ = pi*u0*(diameter)^2/4.0;
   %Fischer et al. 1979, p.345 Rho0=rho fresh water}
   %Lee H.W. and Neville-Jones, 1987, J.Hydr.Eng.,113, 5, p. 616 Rho0=rho_a=rho sea-water}
   gprime = GG*delRho0/Rho0;
   fprintf(1,'delRho0: %5.2f  Rho0: %7.2f gprime: %7.4f',delRho0,Rho0,gprime);
   if abs(dgamdz) > 1e-10 % then
      epsprime = -dgamdz/Rho0;
      % Fischer et al., 1979, page 345 }
      N2 =  GG*epsprime;
      fprintf(1,'N2: %10.6f',N2);        
      ZZstra = konst1*exp(0.25*log(gprime*QQ))...
      /exp(3.0*log(N2)/8.0);
      SSstra = konst2*exp((log(gprime) + ...
      5.*log(ZZstra) - 2.*log(QQ))/3.0);
   else
     ZZstra = DepthDif;
     SSstra = konst3*exp((log(gprime) + ...
     5.*log(ZZstra) - 2.*log(QQ))/3.0);
   end

end% function ANALYTIC }

function [v, derv]=ppvalder(pp,xx)
%   made from ppval by Vlado Malacic (MV) on 2nd October 2019. Besides the values v 
%   at the entries of XX it calculates also the derivatives derv at XX by using 
%   a piecewise polynomial f.  
%   Modifications are in lines with a comment that includes 'MV'.
%   PPVAL  Evaluate piecewise polynomial.
%   V = PPVAL(PP,XX) returns the value, at the entries of XX, of the 
%   piecewise polynomial f contained in PP, as constructed by PCHIP, SPLINE,
%   INTERP1, or the spline utility MKPP.
%
%   V is obtained by replacing each entry of XX by the value of f there.
%   If f is scalar-valued, then V is of the same size as XX. XX may be ND.
%
%   If PP was constructed by PCHIP, SPLINE or MKPP using the orientation of
%   non-scalar function values specified for those functions, then:
%
%   If f is [D1,..,Dr]-valued, and XX is a vector of length N, then V has
%   size [D1,...,Dr, N], with V(:,...,:,J) the value of f at XX(J).
%   If f is [D1,..,Dr]-valued, and XX has size [N1,...,Ns], then V has size
%   [D1,...,Dr, N1,...,Ns], with V(:,...,:, J1,...,Js) the value of f at
%   XX(J1,...,Js).
%   the same holds for derv
%
%   If PP was constructed by INTERP1 using the orientation of non-scalar
%   function values specified for that function, then:
%
%   If f is [D1,..,Dr]-valued, and XX is a scalar, then V has size
%   [D1,...,Dr], with V the value of f at XX.
%   If f is [D1,..,Dr]-valued, and XX is a vector of length N, then V has
%   size [N,D1,...,Dr], with V(J,:,...,:) the value of f at XX(J).
%   If f is [D1,..,Dr]-valued, and XX has size [N1,...,Ns], then V has size
%   [N1,...,Ns,D1,...,Dr], with V(J1,...,Js,:,...,:) the value of f at
%   XX(J1,...,Js).
%
%   Example:
%   Compare the results of integrating the function cos and its spline
%   interpolant:
%
%     a = 0; b = 10;
%     int1 = integral(@cos,a,b);
%     x = a:b; y = cos(x); pp = spline(x,y); 
%     int2 = integral(@(x)ppval(pp,x),a,b);
%
%   int1 provides the integral of the cosine function over the interval [a,b]
%   while int2 provides the integral, over the same interval, of the piecewise
%   polynomial pp that approximates the cosine function by interpolating the
%   computed x,y values.
%
%   Class support for input X and the fields of PP:
%      float: double, single
%
%   See also SPLINE, PCHIP, INTERP1, MKPP, UNMKPP.

%   Carl de Boor
%   Copyright 1984-2015 The MathWorks, Inc.

if isstruct(xx) % we assume that ppval(xx,pp) was used
   temp = xx; xx = pp; pp = temp;
end

%  obtain the row vector xs equivalent to XX
sizexx = size(xx); lx = numel(xx); xs = reshape(xx,1,lx);
%  if XX is row vector, suppress its first dimension
if length(sizexx)==2&&sizexx(1)==1, sizexx(1) = []; end

% take apart PP
[b,c,l,k,dd]=unmkpp(pp);

% for each evaluation site, compute its breakpoint interval
% (mindful of the possibility that xx might be empty)
if lx, [~,index] = histc(xs,[-inf,b(2:l),inf]);
else index = ones(1,lx);
end

% adjust for troubles, like evaluation sites that are NaN or +-inf
infxs = find(xs==inf); if ~isempty(infxs), index(infxs) = l; end
nogoodxs = find(index==0);
if ~isempty(nogoodxs), xs(nogoodxs) = NaN; index(nogoodxs) = 1; end

% now go to local coordinates ...
xs = xs-b(index);

d = prod(dd);
if d>1 % ... replicate xs and index in case PP is vector-valued ...
   xs = reshape(xs(ones(d,1),:),1,d*lx);
   index = d*index; temp = (-d:-1).';
   index = reshape(1+index(ones(d,1),:)+temp(:,ones(1,lx)), d*lx, 1 );
else
   if length(sizexx)>1, dd = []; else dd = 1; end
end

% ... and apply nested multiplication:
v = c(index,1);
derv = (k-1).*c(index,1); % cor MV
% for i=2:k % corr MV 
%    v = xs(:).*v + c(index,i);
% end
for i=2:k-1 % cor MV: for loop only to k-1
   v = xs(:).*v + c(index,i);
   derv = xs(:).*derv + c(index,i).*(k-i); % cor MV
end
i = k; % cor MV - the last iteration for the interpolated value 
v = xs(:).*v + c(index,i); % cor MV


% If evaluating a piecewise constant with more than one piece at NaN, return
% NaN.  With one piece return the constant.
if ~isempty(nogoodxs) && k==1 && l>1
   v = reshape(v,d,lx); v(:,nogoodxs) = NaN;
   derv = reshape(derv,d,lx); derv(:,nogoodxs) = NaN;
end
v = reshape(v,[dd,sizexx]);
derv = reshape(derv,[dd,sizexx]);

if isfield(pp,'orient') && strcmp(pp.orient,'first')
    % spline orientation returns    size(yi) == [d1 ... dk m1 ... mj]
    % but the interp1 usage prefers size(yi) == [m1 ... mj d1 ... dk]
    if ~(isempty(dd) || (isscalar(dd) && dd == 1))
        % The function is non-scalar valued
        if isvector(xx)&&~isscalar(xx)
            permVec = [ndims(v) 1:(ndims(v)-1)];
        else
            ndimsxx = ndims(xx);
            permVec = [(ndims(v)-ndimsxx+1) : ndims(v) 1:(ndims(v)-ndimsxx)];
        end
        v = permute(v,permVec);
        derv = permute(derv,permVec); % cor MV: derv and v are of the same dimension
    end
end 
end % splinder

