How can I calculate this triple integral when I have a rotated euler matrix. All that is attached in the file,

function C_ij = VoigtModelFunction

Cs_kl = [250, 88.3, 87.5, 0, 0, 0 ; ...
    88.3, 250, 87.5, 0, 0, 0 ; ...
    87.5, 87.5, 92, 0, 0, 0 ; ...   %for example(6x6)
    0, 0, 0, 22.7, 0, 0 ; ...
    0, 0, 0, 0, 22.7 0 ; ...
    0, 0, 0, 0, 0, 82.20];
alpha = 0.15235;
if nargin == 0 
    euler = randBunges;
end

euler = euler/180*pi; 
phi1 = euler(1);
phi2  = euler(2);
Phi = euler(3);

alpha=0.15235

Crot= MS_rotEuler(Cs_kl,phi1,Phi,phi2);

VoigtModel = @(phi1, phi2, phi)Crot*ODF_W_function(phi2, alpha).*sin(phi2(i));

C_ij =Crot.*integral3(VoigtModel,0,2*pi,0,2*pi,0,pi);

end
>> VoigtModelFunction

alpha =

      0.1524

CC =

    250.0529   88.2460   87.5006    0.0157   -0.2066   -0.2595
     88.2460  249.9560   87.5318    2.3906    0.0042    0.2677
     87.5006   87.5318   92.0344    0.8356   -0.0703   -0.0026
      0.0157    2.3906    0.8356   22.7319   -0.0006   -0.0967
     -0.2066    0.0042   -0.0703   -0.0006   22.7250    1.2140
     -0.2595    0.2677   -0.0026   -0.0967    1.2140   82.1214

Assignment has more non-singleton rhs dimensions than non-singleton subscripts

Error in ODF_W_function (line 32)
    W(:,ii) = (1./(alpha(ii)*((2*pi).^0.5))) * ...

Error in VoigtModelFunction>@(phi1,phi2,phi)Crot*ODF_W_function(phi2,alpha).*sin(phi2(i)) (line 25)
VoigtModel = @(phi1, phi2, phi)Crot*ODF_W_function(phi2, alpha).*sin(phi2(i));

Error in integral3>@(y,z)FUN(x(1)*ones(size(z)),y,z) (line 138)
            @(y,z)FUN(x(1)*ones(size(z)),y,z), ...

Error in integral2Calc>integral2t/tensor (line 228)
        Z = FUN(X,Y);  NFE = NFE + 1;

Error in integral2Calc>integral2t (line 55)
[Qsub,esub] = tensor(thetaL,thetaR,phiB,phiT);

Error in integral2Calc (line 9)
    [q,errbnd] = integral2t(fun,xmin,xmax,ymin,ymax,optionstruct);

Error in integral3/innerintegral (line 137)
        Q1 = integral2Calc( ...

Error in integralCalc/iterateScalarValued (line 314)
                fx = FUN(t);

Error in integralCalc/vadapt (line 132)
            [q,errbnd] = iterateScalarValued(u,tinterval,pathlen);

Error in integralCalc (line 75)
        [q,errbnd] = vadapt(@AtoBInvTransform,interval);

Error in integral3 (line 121)
    Q = integralCalc(@innerintegral,xmin,xmax,integralOptions);

Error in VoigtModelFunction (line 27)
C_ij =Crot.*integral3(VoigtModel,0,2*pi,0,2*pi,0,pi);

Answer 1

[Mike Hosea]

INTEGRAL3 does not support integrating a matrix of functions. For example, suppose you wanted to integrate

f = @(x,y,z)[cos(x).*sin(y).*z, cos(2*x).*sin(2*y).*2.*z; cos(3*x).*sin(3*y).*3.*z, cos(4*x).*sin(4*y).*4.*z];

over the region x = 0 to pi, y = 0 to pi, and z = 0 to pi. INTEGRAL3 does not do that. INTEGRAL3 requires that if you give your function f inputs of a given size, the output will be the same size. But the function above produces a 2-by-2 matrix when you give it a scalar:

>> f(0,0,0)
ans =
     0     0
     0     0
>> f(pi/4,pi/8,pi/24)
ans =
   0.035421201924390   0.000000000000000
  -0.256543038244258  -0.523598775598299

The only way to use INTEGRAL3 is to split this problem up into the component problems. You can easily integrate f11 = @(x,y,z)cos(x).*sin(y).*z and f12 = @(x,y,z)cos(2*x).*sin(2*y).*2.*z and so forth. Before you attempt to use INTEGRAL3, do some quality control on your integrand function. Make sure it returns the right results, and make sure if the arguments are m-by-n matrices, it returns an m-by-n matrix of results. This must be done element-wise. It cannot matter whether you evaluate the 6 points of a 2-by-3 matrix one at a time by passing in each set of scalar inputs for x, y, and z or all at once by passing in 2-by-3 arrays of inputs for x, y, and z.