# VMT_Rozovskii

## PURPOSE

Computes the frame of reference transpositon as described in Kenworthy

## SYNOPSIS

function [V,log_text] = VMT_Rozovskii(V,A)

## DESCRIPTION

     Computes the frame of reference transpositon as described in Kenworthy
     and Rhoads (1998) ESPL using a Rozovskii type analysis.

     Notes:
         -I extrapolate the velocity at the near surface bin to the water
          surface for the depth averaging (ie, BC at u(z=0) = u(z=bin1))
         -There are cases where the bin corresponding with the bed actually
          contains flow data (i.e., not NaN or zero). For cases where the
          blanking distance DOES exists, I have imposed a BC of U=0 at the bed,
         -In cases where data goes all of the way to the bed, I have divided
          the last bin's velocity by 2, essentially imposing a U=0 at the
          boundary by extrapolating to the bottom of the bin.

     Written by:
     Frank L. Engel, USGS (fengel@usgs.gov)
     Last edited: F.L. Engel, USGS, 2/20/2013

## CROSS-REFERENCE INFORMATION

This function calls:

This function is called by:

- VMT_ProcessTransects Driver program to process multiple transects at a single cross-section

## SOURCE CODE

    0001 function [V,log_text] = VMT_Rozovskii(V,A)
    0002 % Computes the frame of reference transpositon as described in Kenworthy
    0003 % and Rhoads (1998) ESPL using a Rozovskii type analysis.
    0004 %
    0005 % Notes:
    0006 %     -I extrapolate the velocity at the near surface bin to the water
    0007 %      surface for the depth averaging (ie, BC at u(z=0) = u(z=bin1))
    0008 %     -There are cases where the bin corresponding with the bed actually
    0009 %      contains flow data (i.e., not NaN or zero). For cases where the
    0010 %      blanking distance DOES exists, I have imposed a BC of U=0 at the bed,
    0011 %     -In cases where data goes all of the way to the bed, I have divided
    0012 %      the last bin's velocity by 2, essentially imposing a U=0 at the
    0013 %      boundary by extrapolating to the bottom of the bin.
    0014 %
    0015 % Written by:
    0016 % Frank L. Engel, USGS (fengel@usgs.gov)
    0017 % Last edited: F.L. Engel, USGS, 2/20/2013
    0018
    0019 %disp('Performing Rozovskii analysis...')
    0020 log_text = {'      Performing Rozovskii analysis...'};
    0021
    0022 % Calculate dy and dz for each meaurement point
    0023 dy=mean(diff(V.mcsDist(1,:)));%m
    0024 dz=mean(diff(V.mcsDepth(:,1)));%m
    0025
    0026 % Pull the velocity data, and convert NaNs to zeros (needed for proper
    0027 % depth integration to the boundary)
    0028 u = V.u; v = V.v;
    0029 idx = isnan(u) | isnan(v);
    0030 u(idx) = 0; v(idx) = 0;
    0031
    0032 % Pull depth data
    0033 d = V.mcsDepth;
    0034
    0035 % Pull the bed data
    0036 b = V.mcsBed;
    0037
    0038 % Work on each vertical (not necesarily the same as ensemble) seperately
    0039 for i = 1:size(u,2)
    0040     % Finds closest cell to bed
    0041     [~, array_position(i)]...
    0042         = min(abs(d(:,i) - b(i)));
    0043     for j = 1:array_position(i)
    0044         if j == 1 % Near water surface
    0045             % Compute first bin by exprapolating velocity to the water
    0046             % surface. WSE = 0. Imposes BC u(z=0) = u(z=bin1)
    0047             du_i(j,i) = u(j,i)*(d(j+1,i)-d(j,i))...
    0048                 + u(j,i)*(d(j,i)-dz/2-0);
    0049             dv_i(j,i) = v(j,i)*(d(j+1,i)-d(j,i))...
    0050                 + v(j,i)*(d(j,i)-dz/2-0);
    0051         elseif j < array_position(i) % Inbetween
    0052             du_i(j,i) = u(j,i)*(d(j+1,i)-d(j-1,i))/2;
    0053             dv_i(j,i) = v(j,i)*(d(j+1,i)-d(j-1,i))/2;
    0054         elseif j == array_position(i) % Near bed
    0055             indx = find(u(:,i) ~= 0);  %Revision PRJ 9-1-09
    0056             if isempty(indx)
    0057                 du_i(:,i) = NaN;
    0058                 dv_i(:,i) = NaN;
    0059             else
    0060                 l = indx(end);
    0061                 k = j - l;
    0062                 % Computes du from last good bin to the bed by linear
    0063                 % interpolation. IMPOSES BC: u=0 at the bed
    0064                 % Paints everything below last bin as NaN
    0065                 du_i(j-k+2:size(u,2),i) = NaN;
    0066                 dv_i(j-k+2:size(u,2),i) = NaN;
    0067             end
    0068         end
    0069     end
    0070
    0071     % Depth averaged quantities
    0072     U(i) = nansum(du_i(:,i))/d(array_position(i),i);
    0073     V1(i) = nansum(dv_i(:,i))/d(array_position(i),i);
    0074     U_mag(i) = sqrt(U(i)^2+V1(i)^2); % resultant vector
    0075
    0076     % Angle of resultant vector from a perpendicular line along the
    0077     % transect
    0078     phi(i) = atan(V1(i)/U(i));
    0079     phi_deg(i) = phi(i).*180/pi;
    0080
    0081     % Compute Rozovskii variables at each bin
    0082     for j = 1:array_position(i)
    0083         uu(j,i) = sqrt(u(j,i)^2+v(j,i)^2);
    0084         if (u(j,i) < 0) && (v(j,i) < 0)
    0085             theta(j,i) = atan(v(j,i)/u(j,i)) - pi();
    0086         elseif (u(j,i) < 0) && (v(j,i) > 0)
    0087             theta(j,i) = atan(v(j,i)/u(j,i)) + pi();
    0088         else
    0089             theta(j,i) = atan(v(j,i)/u(j,i));
    0090         end
    0091         theta_deg(j,i) = theta(j,i).*180/pi;
    0092         up(j,i) = uu(j,i)*cos(theta(j,i)-phi(i));
    0093         us(j,i) = uu(j,i)*sin(theta(j,i)-phi(i));
    0094         upy(j,i) = up(j,i)*sin(phi(i));
    0095         upx(j,i) = up(j,i)*cos(phi(i));
    0096         usy(j,i) = us(j,i)*cos(phi(i));
    0097         usx(j,i) = us(j,i)*sin(phi(i));
    0098         depths(j,i) = d(j,i);
    0099
    0100         % Compute d_us to check for zero secondary discharge constraint
    0101         if j == 1 % Near water surface
    0102             dus_i(j,i) = us(j,i)*(d(j+1,i)-d(j,i))...
    0103                 + us(j,i)*(d(j,i)-dz/2-0);
    0104         elseif j < array_position(i) % Inbetween
    0105             dus_i(j,i) = us(j,i)*(d(j+1,i)-d(j-1,i))/2;
    0106         end
    0107         % Sum dus_i - this should be near zero for each ensemble
    0108         q_us(i) = nansum(dus_i(:,i));
    0109     end
    0110
    0111     % Resize variables to be the same as input grids
    0112     uu(j+1:size(u,1),i) = NaN;
    0113     theta(j+1:size(u,1),i) = NaN;
    0114     theta_deg(j+1:size(u,1),i) = NaN;
    0115     up(j+1:size(u,1),i) = NaN;
    0116     us(j+1:size(u,1),i) = NaN;
    0117     upy(j+1:size(u,1),i) = NaN;
    0118     usy(j+1:size(u,1),i) = NaN;
    0119     upx(j+1:size(u,1),i) = NaN;
    0120     usx(j+1:size(u,1),i) = NaN;
    0121     depths(j+1:size(u,1),i) = NaN;
    0122     dus_i(j+1:size(u,1),i) = NaN;
    0123 end
    0124
    0125 % Display error message if rozovskii computation of q_us doesn't sum to
    0126 % zero
    0127 if q_us > 1e-4
    0128     %disp('Warning: Rozovskii secondary velocities not satisfying continuity!')
    0129     log_text = vertcat(log_text,...
    0130         '      Warning: Rozovskii secondary velocities',...
    0131         '      not satisfying continuity!');
    0132 else
    0133     %disp('Computation successfull: Rozovskii secondary velocities satisfy continuity')
    0134     log_text = vertcat(log_text,...
    0135         '      Computation successfull: Rozovskii',...
    0136         '      secondary velocities',...
    0137         '      satisfy continuity.');
    0138 end
    0139
    0140 % Rotate local velocity vectors into global coordinate system by
    0141 % determining the angle of the transect using endpoint locations. The
    0142 % function "vrotation" is a standard rotation matrix
    0143 XStheta = atan((V.mcsY(1,end)-V.mcsY(1,1))/(V.mcsX(1,end)-V.mcsX(1,1)));
    0144 XSalpha = XStheta - pi/2;
    0145 [ux, uy, uz] = vrotation(V.u,V.v,V.w,XSalpha);
    0146
    0147
    0148 % Write outputs to Roz substructure
    0149 V.Roz.U = U;
    0150 V.Roz.V = V1;
    0151 V.Roz.U_mag  = U_mag;
    0152 V.Roz.phi = phi;
    0153 V.Roz.phi_deg = phi_deg;
    0154 V.Roz.u = V.u;
    0155 V.Roz.v = V.v;
    0156 V.Roz.u_mag = u;
    0157 V.Roz.depth = depths;
    0158 V.Roz.theta = theta;
    0159 V.Roz.theta_deg = theta_deg;
    0160 V.Roz.up = up;
    0161 V.Roz.us = us;
    0162 V.Roz.upy = upy;
    0163 V.Roz.usy = usy;
    0164 V.Roz.upx = upx;
    0165 V.Roz.usx = usx;
    0166 V.Roz.ux = ux;
    0167 V.Roz.uy = uy;
    0168 V.Roz.uz = uz;
    0169 V.Roz.alpha = XSalpha;
    0170
    0171 %disp('Rozovskii analysis complete. Added .Roz structure to V data structure.')
    0172 log_text = vertcat(log_text, '      Rozovskii analysis complete.');
    0173 % directory = pwd;
    0174 % fileloc = [directory '\' filename '.mat'];
    0175 % disp(fileloc)