//JLJ
//Sweep space for speed
//Flowrate calculation and visualitation
clear;
D=0.1;// [m] Internal diameter of pipe
R=D/2;//[m] Internal radius of pipe
in=100; //Grid spacing, xcord
jn=100; //Grid spacing, ycord
vmax=2;// [m/s] Max flow speed, centerline speed (x=y=0)
dA=D/in*D/jn;//[m^2]
//Physical multiples point measurement data
ns=5; //Number of different measurements, equidistant stations
xp=linspace(-R,R,in);
yp=linspace(-R,R,jn);
flowrateideal=0;
flowratemultipoint=0;
dadda=0;
i=1;
j=1;
for i=1:in
    for j=1:jn
       r=((xp(i)^2+yp(j)^2)^0.5)
       if r>R then
           r=R;
        end
      //  theta=atan2(yp(j),xp(i))
        v(i,j)=vmax*(1-r/R)^(1/7) //Speed magnitude
        if r<R then
        flowrateideal=flowrateideal+v(i,j)*dA //Flowrate calculation
        dadda=dadda+dA
        end
    end
end

unohandler=figure(1,"BackgroundColor",[1 1 1])
plot3d1(xp,yp,v)
zlabel("Velocity Magnitude m/s")
xtitle("Turbolent velocity profile in a "+sprintf("%.4f",R)+ " m radius pipe, vmax "+ sprintf("%.1f",vmax)+" m/s Flowrate "+sprintf("%f",flowrateideal)+" m^3/s" )
//Calculation of physical measurement
duehandler=figure(2,"BackgroundColor",[1 1 1])
ip=1;
for ip=1:ns
       r=R/ns*ip-R/ns*0.5
   //    if ip==ns then
   //    r=R/ns*ip-R/ns*0.5
   //    end
       ru=R/ns*ip
       rd=R/ns*(ip-1)
       if ip==1 then
       rd=0;
       end
       vmp(i)=vmax*(1-r/R)^(1/7) //Speed magnitude
       flowratemultipoint=flowratemultipoint+vmp(i)*(3.14*(ru^2-rd^2)) //Flowrate calculation
       for ii=1:in
           for jj=1:jn
               if ((((xp(ii)^2+yp(jj)^2))^0.5)<r) then
               vout(ii,jj)=vmp(i);
               else
               vout(ii,jj)=0
               end
            end
       end
      plot3d1(xp,yp,vout)
      zlabel("Velocity Magnitude m/s")
      xtitle("Turbolent velocity profile in a "+sprintf("%.4f",R)+ " m radius pipe, vmax "+ sprintf("%.1f",vmax)+" m/s Flowrate "+sprintf("%f",flowratemultipoint)+" m^3/s" )
end


