Organizing and double integrating a large function

I am trying to build a function using nested functions rather than trying to write out the whole function (sure to make errors). Then I would like to integrate that function with respect to r and psi. How do I set this up so I can have visability to the functions inside while still being able to integrate it.

function [dH] = dHf(r,psi)
% set constants
rho = 1.225;
Xin = 0.3048;
Yin = 0;

R = 0.1524;

[Cl, Cd, Cm, Re_values] = nacaData('4412');

% Velocities at the hub
% p = roll rate, q = pitch rate, r = yaw rate
u = 0 ;v = 0;w = 0;p = 0;q = 0;yawr = 0;
% North Rotor
Vhubxn = u + Yin*yawr;
Vhubyn = v + Xin*yawr;
Vhubzn = w - Yin*p - Xin*q;

Omega = 1100;
lambda = 0.5;

    % Define dH based on dFx 
    dFx = @Fx;
    dH = (-1).*dFx(r,psi).*sin(psi);
        
    % Define dFx based on phii, dL and dD 
        function [dFx] = Fx(r,psi)
        dD = @dDr;
        dL = @dLi;
        phii = @phi;
        dFx = dL(r,psi)*sin(phii(r,psi)) + dD(r,psi)*cos(phii(r,psi));
        % define phii   
            function [Phi] = phi(r,psi)
            UP = @up;
            UT = @ut;
            Phi = atan(UP(r,psi)/UT(r,psi));
                function [upo] = up(r,psi)
                upo = Omega*R* lambda + Vhubzn - ((-1).^0) * r * rho * sin(psi) - r*q*cos(psi);
                end
                function [uto] = ut(r,psi)
                uto = Omega*r + Vhubxn*sin(psi) + ((-1).^0) * Vhubyn *cos(psi);
                end
        
            end
          % Drag definition ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
            function [D] = dDr(r,psi)
            UP = @up;
            UT = @ut;
            c = @cb;
            Cdr = @cDo;
            
            
            D = 0.5*Cdr(r,psi)*rho*(UT(r,psi).^2+UP(r,psi).^2)*c(r);
                
                function [CD] = cDo(r,psi)
                alpha = @a;
                Cda = Cd(95:201,3);
                Cdaa = Cd(95:201,2);
                Cdeq = polyfit(Cdaa,Cda,3);

                CD = Cdeq(1,1)*((alpha(r,psi)).^3) + Cdeq(1,2)*((alpha(r,psi)).^2) + Cdeq(1,3)*(alpha(r,psi)) + Cdeq(1,4);
                    
                    function [ao] = a(r,psi)
                    thetab = @theta;
                    phi = @phii;    
                    
                    ao = thetab(r) - phi(r,psi);
                        function [thetablade] = theta(r)
                        thetablade = -1.1896*r+0.375;
                        end

                        function [Phi] = phii(r,psi)
                        UP = @up;
                        UT = @ut;
                        Phi = atan(UP(r,psi)/UT(r,psi));
                            function [upo] = up(r,psi)
                            upo = Omega*R* lambda + Vhubzn - ((-1).^0) * r * rho * sin(psi) - r*q*cos(psi);
                            end
                            function [uto] = ut(r,psi)
                            uto = Omega*r + Vhubxn*sin(psi) + ((-1).^0) * Vhubyn *cos(psi);
                            end
        
                        end


                    end

                end
                
                function [upo] = up(r,psi)
                upo = Omega*R* lambda + Vhubzn - ((-1).^0) * r * rho * sin(psi) - r*q*cos(psi);
                end
                function [uto] = ut(r,psi)
                uto = Omega*r + Vhubxn*sin(psi) + ((-1).^0) * Vhubyn *cos(psi);
                end
                
                function [C] = cb(r)
                C = -0.12467*r+0.031;
                end

            end
%~~ Lift ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~      
            function [L] = dLi(r,psi)
            UP = @up;
            UT = @ut;
            c = @cb;
            Clo = @cLo;
            
            L = 0.5*Clo(r,psi)*rho*(UT(r,psi).^2+UP(r,psi).^2)*c(r);   
%~~~~~ Cl ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~               
                function [CL] = cLo(r,psi)
                alpha = @a;
                Cla = Cl(95:201,3);
                Claa = Cl(95:201,2);
                Cleq = polyfit(Claa,Cla,3);

                CL = Cleq(1,1)*((alpha(r,psi)).^3) + Cleq(1,2)*((alpha(r,psi)).^2) + Cleq(1,3)*(alpha(r,psi)) + Cleq(1,4);
                    
                    function [ao] = a(r,psi)
                    thetab = @theta;
                    phi = @phii;    
                    
                    ao = thetab(r) - phi(r,psi);
                        function [thetablade] = theta(r)
                        thetablade = -1.1896*r+0.375;
                        end

                        function [Phi] = phii(r,psi)
                        UP = @up;
                        UT = @ut;
                        Phi = atan(UP(r,psi)/UT(r,psi));
                            function [upo] = up(r,psi)
                            upo = Omega*R* lambda + Vhubzn - ((-1).^0) * r * rho * sin(psi) - r*q*cos(psi);
                            end
                            function [uto] = ut(r,psi)
                            uto = Omega*r + Vhubxn*sin(psi) + ((-1).^0) * Vhubyn *cos(psi);
                            end
        
                        end


                    end 
                end
%~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
                function [upo] = up(r,psi)
                upo = Omega*R* lambda + Vhubzn - ((-1).^0) * r * rho * sin(psi) - r*q*cos(psi);
                end
                function [uto] = ut(r,psi)
                uto = Omega*r + Vhubxn*sin(psi) + ((-1).^0) * Vhubyn *cos(psi);
                end
             
                function [C] = cb(r)
                C = -0.12467*r+0.031;
                end
            
            
            
            end 
    


          end
end

[Torsten]

How do I set this up so I can have visability to the functions inside while still being able to integrate it.

I don't understand what you mean. Program the function you are trying to integrate in a way so that you can best follow the written code.

If you want to use "integral2" for integration, your function should accept matrices for r and psi and return f(r,psi) as a matrix of the same size.

[Jeff Coulter]

When I try to use integral2 or trapz, it doesn't see what the internal equations are that build the top level function. it returns not enough input arguments because it doesn't see past the function

example of what I am thinking;

integral2(A, x1,x2,y1,y2)

A = B*C;

C = x *D;

D = y * 2x;

B = y*x;

[Torsten]

If you like such a structure, use it. I don't.

integral2(@A, x1,x2,y1,y2)

function valueA = A(x,y)
  valueA = B(x,y).*C(x,y)
end
function valueB = B(x,y)
  valueB = y.*x;
end
function valueC = C(x,y)
  valueC = x.*D(x,y);
end
function valueD = D(x,y)
  valueD = y.*2x;
end