source: issm/trunk-jpl/examples/EsaWahr/wahr.m@ 22793

% compute semi-analytic solutions for a disc load 
function [vert, horz] = wahr(disc_rad,xi,love_h,love_l); 

        disc_rad = disc_rad/1000; % km 
        % compute P(x), dP(x)/dx, d2P(x)/dx2
        %---------------------------------------------------------------------
        % compute p_value 
        theta=km2deg(xi/1000)';
        ang = theta/180*pi; 
        alpha=cos(ang);
        m=length(alpha);
        n=length(love_h)-1; 
        p_value = p_polynomial_value(m,n,alpha);
        p_prime = p_polynomial_prime(m,n,alpha);
        %---------------------------------------------------------------------
        nn=[0:n];
        nn_plus_1=nn+1; 

        % disc radius in degree 
        disc_rad = km2deg(disc_rad)/180*pi; 
        tau=zeros(size(love_h)); 
        tau(1) = 0.5*(1-cos(disc_rad)); % tau_0 
        p_value_disc = p_polynomial_value(1,n+1,cos(disc_rad));
        p_prime_disc = p_polynomial_prime(1,n,cos(disc_rad));
        for jj=2:n+1
                nnn = jj-1; 
                tau(jj) = 0.5 * (p_value_disc(jj-1) - p_value_disc(jj+1)); 
        end

        const=zeros(size(love_h)); 
        for jj=1:n+1
                const(jj) = 1/(2*(jj-1)+1); 
        end

        disc=sum(bsxfun(@times,p_value,tau'),2); 

        g1 = -sum(bsxfun(@times,p_value,(tau.*love_h.*const)'),2); 
        g5 = -sum(bsxfun(@times,-sin(ang),bsxfun(@times,p_prime,(tau.*love_l.*const)')),2); 

        % coeff 
        coeff = 1000*4*pi*(6.67408*10^-11)*(6.3781*10^6)/9.81; 

        % vertical and horizontal solutions in mm 
        vert = g1*coeff*1000; % mm 
        horz = g5*coeff*1000; % mm