title: speed up code - vectorization
body: <p>I'm not really familiar with <code>vectorization</code>, but I am aware that, amongst <code>MATLAB</code>'s strengths, code vectorization is probably the most rewarded.</p>
  
  <p>I have this code:</p>
  
  <pre><code>ikx= (-Nx/2:Nx/2-1)*dk1;
  iky= (-Ny/2:Ny/2-1)*dk2;
  ikz= (-Nz/2:Nz/2-1)*dk3;
  
  [k1,k2,k3] = ndgrid(ikx,iky,ikz);
  k = sqrt(k1.^2 + k2.^2 + k3.^2);
  Cij = zeros(3,3,Nx,Ny,Nz);
  count = 0;
  for ii = 1:Nx
      for jj = 1:Ny
          for kk = 1:Nz
              if ~isequal(k1(ii,jj,kk),0)
                  count = count +1;
                  fprintf('iteration step %i\r\n',count)
                  E_int = interp1(k_vec,E_vec,k(ii,jj,kk),'spline','extrap');
                  beta = c*gamma./(k(ii,jj,kk).*sqrt(E_int));
                  k30 = k3(ii,jj,kk) + beta*k1(ii,jj,kk);
                  k0 = sqrt(k1(ii,jj,kk)^2 + k2(ii,jj,kk)^2 + k30^2);
                  Ek0 = 1.453*(k0^4/((1 + k0^2)^(17/6)));
                  B = sigmaiso*sqrt((Ek0./(k0.^2))*((dk1*dk2*dk3)/(4*pi)));
                  C1 = ((beta.*k1(ii,jj,kk).^2).*(k0.^2 - 2*k30.^2 + k30.*beta.*k1(ii,jj,kk)))./(k(ii,jj,kk).^2.*(k1(ii,jj,kk).^2 + k2(ii,jj,kk).^2));
                  C2 = ((k2(ii,jj,kk).*(k0.^2))./((k1(ii,jj,kk).^2 + k2(ii,jj,kk).^2).^(3/2))).*atan2((beta.*k1(ii,jj,kk).*sqrt(k1(ii,jj,kk).^2 + k2(ii,jj,kk).^2)),(k0.^2 - k30.*beta.*k1(ii,jj,kk)));
                  xhsi1 = C1 - C2.*(k2(ii,jj,kk)./k1(ii,jj,kk));
                  xhsi2 = C1.*(k2(ii,jj,kk)./k1(ii,jj,kk)) + C2;
                  Cij(1,1,ii,jj,kk) = B.*((k2(ii,jj,kk).*xhsi1)./(k0));
                  Cij(1,2,ii,jj,kk) = B.*((k3(ii,jj,kk)-k1(ii,jj,kk).*xhsi1+beta.*k1(ii,jj,kk))./(k0));
                  Cij(1,3,ii,jj,kk) = B.*(-k2(ii,jj,kk)./(k0));
                  Cij(2,1,ii,jj,kk) = B.*((k2(ii,jj,kk).*xhsi2-k3(ii,jj,kk)-beta.*k1(ii,jj,kk))./(k0));
                  Cij(2,2,ii,jj,kk) = B.*((-k1(ii,jj,kk).*xhsi2)./(k0));
                  Cij(2,3,ii,jj,kk) = B.*(k1(ii,jj,kk)./(k0));
                  Cij(3,1,ii,jj,kk) = B.*(k2(ii,jj,kk).*k0./(k(ii,jj,kk).^2));
                  Cij(3,2,ii,jj,kk) = B.*(-k1(ii,jj,kk).*k0./(k(ii,jj,kk).^2));               
              end
          end
      end
  end
  </code></pre>
  
  <p>Generally, I might avoid the nested <code>for</code> loops; nonetheless, the <code>if</code> statement on <code>k1</code> values is currently directing me towards the classical and old-fashion code structure.</p>
  
  <p>I blatantly would like to bypass the presence of the <code>for</code> loops in favour of vectorized and more elegant solution.</p>
  
  <p>Any support is more than welcome.</p>
  
  <p><strong>EDIT</strong></p>
  
  <p>To let better understand what the code is expected to perform, I hereby provide you with some basics:</p>
  
  <p><img src="https://i.stack.imgur.com/Olowo.png" alt="enter image description here"></p>
  
  <p><img src="https://i.stack.imgur.com/yB2fZ.png" alt="enter image description here"></p>
  
  <p><img src="https://i.stack.imgur.com/RGB4J.png" alt="enter image description here"></p>
  
  <p><img src="https://i.stack.imgur.com/bMZe4.png" alt="enter image description here"></p>
  
  <p><strong>EDIT2</strong></p>
  
  <p>As @Floris advised, I came up with this alternative solution:</p>
  
  <pre><code>ikx= (-Nx/2:Nx/2-1)*dk1;
  iky= (-Ny/2:Ny/2-1)*dk2;
  ikz= (-Nz/2:Nz/2-1)*dk3;
  
  [k1,k2,k3] = ndgrid(ikx,iky,ikz);
  k = sqrt(k1.^2 + k2.^2 + k3.^2);
  
  ii = (ikx ~= 0);
  k1w = k1(ii,:,:);
  k2w = k2(ii,:,:);
  k3w = k3(ii,:,:);
  kw = k(ii,:,:);
  
  E_int = interp1(k_vec,E_vec,kw,'spline','extrap');
  beta = c*gamma./(kw.*sqrt(E_int));
  
  k30 = k3w + beta.*k1w;
  k0 = sqrt(k1w.^2 + k2w.^2 + k30.^2);
  Ek0 = (1.453*k0.^4)./((1 + k0.^2).^(17/6));
  B = sqrt((2*(pi^2)*(l^3))*(Ek0./(V*k0.^4)));
  
  k1w_2 = k1w.^2;
  k2w_2 = k2w.^2;
  k30_2 = k30.^2;
  k0_2 = k0.^2;
  kw_2 = kw.^2;
  
  C1 = ((beta.*k1w_2).*(k0_2 - 2.*k30_2 + beta.*k1w.*k30))./(kw_2.*(k1w_2 + k2w_2));
  C2 = ((k2w.*k0_2)./((k1w_2 + k2w_2).^(3/2))).*atan2((beta.*k1w).*sqrt(k1w_2 + k2w_2),(k0_2 - k30.*k1w.*beta));
  
  xhsi1 = C1 - (k2w./k1w).*C2;
  xhsi2 = (k2w./k1w).*C1 + C2;
  
  Cij = zeros(3,3,Nx,Ny,Nz);
  
  Cij(1,1,ii,:,:) = B.*(k2w.*xhsi1);
  Cij(1,2,ii,:,:) = B.*(k3w - k1w.*xhsi1 + beta.*k1w);
  Cij(1,3,ii,:,:) = B.*(-k2w);
  Cij(2,1,ii,:,:) = B.*(k2w.*xhsi2 - k3w - beta.*k1w);
  Cij(2,2,ii,:,:) = B.*(-k1w.*xhsi2);
  Cij(2,3,ii,:,:) = B.*(k1w);
  Cij(3,1,ii,:,:) = B.*((k0_2./kw_2).*k2w);
  Cij(3,2,ii,:,:) = B.*(-(k0_2./kw_2).*k1w);
  </code></pre>