function [h_o,index_1]=orth_func(nx,ny)
%
% [h_o,index_1]=orth_func(nx,ny)
%
% function to be used to compute 20 orthonormal functions which are the
% important part to compute the image spectrum parameters
% input: [nx,ny] are the rows and cols of image spectrum
% output: 
%    h_o is the orthnonrmal functions 
%    index_1 is the flag showed the area between the maximum wavenumber and
%    minimum wave number. 
%

    dkx=2.0*pi./(nx*16.0);                        % the wavenumber spacing of range direction
    dky=2.0*pi./(ny*20.0);                        % the wavenumber spacing of azimuth direction

    kmin=2.0*pi/624;                              % the minimum wavenumber
    kmax=2.0*pi/60;                               % the maximum wavenumber

    gamma=2.0;                                    % azimuth bunching parameter

    einsx=ones(nx,1);
    einsy=ones(ny,1);

    % build the wavenumber matrix which shows the wavenumber of every pixel
    akx=dkx.*((1:nx)'*einsy'-fix(nx/2.0));
    aky=(dky.*(einsx*(1:ny)-fix(ny/2.0))')';
    akmod=sqrt(akx.^2+aky.^2);

    dummy_a=(gamma^2-gamma^4)/(gamma^2*kmin^2-kmax^2);              % the parameter : a1
    dummy_b=(kmax^2-gamma^4*kmin^2)/(kmax^2-gamma^2*kmin^2);        % the parameter : a2
    dummy_c=(dummy_a*akx.^4+dummy_b*akx.^2+aky.^2);

    alpha2=atan(aky./akx);                                          % the paramater : alpha_phi
    alpha2(nx/2,ny/2)=alpha2(nx/2-1,ny/2);
    alpha1=2.0*(log(sqrt(dummy_c))-log(kmin))/(log(kmax)-log(kmin))-1;% the parameter : alpha_k

    g_o=zeros(nx,ny,4);
    f_o=zeros(nx,ny,5);
    
    % set the range between the maximum wavenumber and minimum wavenumber
    index_0 = (alpha1< -1.) | (alpha1 > 1.) | (alpha2 <-pi/2.)|(alpha2>pi/2.);
    index_1 = (alpha1>=-1.) & (alpha1 <=1.) & (alpha2>=-pi/2.)&(alpha2<=pi/2.);
    
    % compute 4 Gegenbauer polynomials
    alpha1(index_0)=nan;
    g_o(:,:,1)=ones(nx,ny)*sqrt(3.0/4.0).*(1-alpha1.^2).^0.5;
    g_o(:,:,2)=sqrt(5.0/12.0)*3*alpha1.*(1-alpha1.^2).^0.5;
    g_o(:,:,3)=sqrt(7.0/24.0)*(15.0/2.0*alpha1.^2-3.0/2.0).*(1-alpha1.^2).^0.5;
    g_o(:,:,4)=sqrt(9.0/40.0)*(35.0/2.*alpha1.^3-15.0/2.*alpha1).*(1-alpha1.^2).^0.5;
    
    % compute 5 harmonic functions
    f_o(:,:,1)=ones(nx,ny)*sqrt(1.0/pi);
    f_o(:,:,2)=sqrt(2./pi)*sin(2.0*alpha2);
    f_o(:,:,3) = sqrt(2./pi)*cos(2.0*alpha2);
    f_o(:,:,4) = sqrt(2./pi)*sin(4.0*alpha2);
    f_o(:,:,5) = sqrt(2./pi)*cos(4.0*alpha2);
    
    % compute the weight function
    weight = sqrt( 2.0*(dummy_b*akx.^2+2.0*dummy_a*akx.^4+aky.^2)./ ...
             ( (akx.^2+aky.^2).*(dummy_b*akx.^2+dummy_a*akx.^4+aky.^2)*(log(kmax)-log(kmin))));
    weight(index_0)=nan;    
    
    % compute 20 orthonormal functions
    h_o=zeros(nx,ny,20);
    for i=1:4
        for j=1:5
            dummy = f_o(:,:,j).*g_o(:,:,i).*weight;
            dummy(index_0) = nan;
            h_o(:,:,(i-1)*5+j) =  dummy;
        end
    end
end