title: Integral Calculation
body: <p>I would like to perform what follows:</p>
  
  <p><img src="https://i.stack.imgur.com/HsFla.jpg" alt="enter image description here"></p>
  
  <p>where</p>
  
  <p><img src="https://i.stack.imgur.com/IfGxH.jpg" alt="enter image description here"></p>
  
  <p>The parameters shown in the latter figure can be obtained as follows:</p>
  
  <pre><code>%% Inizialization
  time = 614.4; % Analysis Time
  Uhub = 11;
  HubHt = 90;
  alpha = 0.14;
  TI = 'A'; % Turbulent Intensity (A,B,C as in the IEC or Specific value)
  N1 = 4096;
  N2 = 32;
  N3 = 32;
  
  N = N1*N2*N3; % Total Number of Point
  t = 0:(time/(N1-1)):time; % Sampled Time Vector
  
  L1 = Uhub*time; % Box length along X
  L2 = 150; % Box length along Y
  L3 = 220; % Box length along Z
  
  dx = L1/N1; % Grid Resolution along X-axis
  dy = L2/N2; % Grid Resolution along Y-axis
  dz = L3/N3; % Grid Resolution along Z-axis
  
  V = L1*L2*L3; % Analysis Box Volume
  
  gamma = 3.9; % Turbulent Eddies Distorsion Factor
  c = 1.476;
  b = 5.6;
  
  if HubHt &lt; 60
      lambda1 = 0.7*HubHt;
  else
      lambda1 = 42;
  end
  L = 0.8*lambda1;
  
  if isequal(TI,'A')
      Iref = 0.16;
      sigma1 = Iref*(0.75*Uhub + b);
  elseif isequal(TI,'B')
      Iref = 0.14;
      sigma1 = Iref*(0.75*Uhub + b);
  elseif isequal(TI,'C')
      Iref = 0.12;
      sigma1 = Iref*(0.75*Uhub + b);
  else
      sigma1 = str2num(TI)*Uhub/100;
  end
  
  sigma_iso = 0.55*sigma1;
  sigma2 = 0.7*sigma1;
  sigma3 = 0.5*sigma1;
  
  
  %% Wave number vectors
  ik1 = cat(2,(-N1/2:-1/2),(1/2:N1/2));
  ik2 = -N2/2:N2/2-1;
  ik3 = -N3/2:N3/2-1;
  
  [x y z] = ndgrid(ik1,ik2,ik3);
  
  k1 = reshape((2*pi*L/L1)*x,N1*N2*N3,1);
  k2 = reshape((2*pi*L/L2)*y,N1*N2*N3,1);
  k3 = reshape((2*pi*L/L3)*z,N1*N2*N3,1);
  
  k = sqrt(k1.^2 + k2.^2 + k3.^2);
  
  %% Calculation of beta by means of the Energy Spectrum Integration
  E = @(k) (1.453*k.^4)./((1 + k.^2).^(17/6));
  
  %//Independent integration on segments 
  Local_int = arrayfun(@(i)quad(E,i-1,i), 2:(N1*N2*N3));      
  %//integral additivity + cumulative removal of queues
  E_int = 1.5 - [0 fliplr(cumsum(fliplr(Local_int)))]; %//To remove queues
  
  E_int = reshape(E_int,N,1);
  S = k.*sqrt(E_int);
  beta = (c*gamma)./S;
  
  %% Help Parameters
  k30 = k3 + beta.*k1;
  
  k0 = sqrt(k1.^2 + k2.^2 + k30.^2);
  
  C1 = (beta.*k1.^2.*(k1.^2 + k2.^2 - k3.*k30))./(k.^2.*(k1.^2 + k2.^2));
  C2 = (k2.*k0.^2./((k1.^2 + k2.^2).^(3/2))).*atan2((beta.*k1.*sqrt(k1.^2 + k2.^2)),(k0.^2 - k30.*k1.*beta));
  
  xhsi1 = C1 - (k2./k1).*C2;
  xhsi2 = (k2./k1).*C1 + C2;
  
  E_k0 = (1.453*k0.^4)./((1 + k0.^2).^(17/6));
  </code></pre>
  
  <p>For instance, typing in</p>
  
  <pre><code>phi_33 = @(k2,k3) (E_k0./(4*pi.*k.^4)).*((k1.^2 + k2.^2));
  
  F_33 = arrayfun(@(i) dblquad(phi_33,k3(i),k3(i+1),k2(i),k2(i+1)), 1:((N1*N2*N3)-1));
  </code></pre>
  
  <p>Matlab retrieves the following error msg:</p>
  
  <pre><code>Error using  + 
  Matrix dimensions must agree.
  
  Error in @(k2,k3)(E_k0./(4*pi.*k.^4)).*((k1.^2+k2.^2))
  </code></pre>
  
  <p>Do you have a clue how to overcome this issue?</p>
  
  <p>I really look forward to hearing from you.</p>
  
  <p>Best regards,
  FPE</p>