Energy Spectrum Analysis from vorticity in the spectrum space

gistfile1.matlab

classdef SpectrumAnalysis < handle
    properties (SetAccess = public, GetAccess = public)
        EperWavenumber
        wavenumber
        energyContainWavenumber
        t_start
        t_end
        p_lowb
        p_upperb
        numFourier
        numSpherical
    end
    
    properties (SetAccess = private, GetAccess = private)
        vors_real
        vors_imag
    end
        
    methods
        function obj = SpectrumAnalysis(filename,...
                startT,endT,p_low,p_upperb)
            ncid = netcdf.open(filename,'NC_NOWRITE');
            varid_real = netcdf.inqVarID(ncid,'vors_real');
            varid_imag = netcdf.inqVarID(ncid,'vors_imag');
            obj.vors_real = netcdf.getVar(ncid,varid_real);
            obj.vors_imag = netcdf.getVar(ncid,varid_imag);
            netcdf.close(ncid)
            
            if nargin == 1
                obj.t_start = 300;
                obj.t_end =size(obj.vors_real,4);
                obj.p_lowb = floor(0.25*size(obj.vors_real,3));
                obj.p_upperb = floor(0.75*size(obj.vors_real,3));
                obj.numFourier = size(obj.vors_real,1);
                obj.numSpherical = size(obj.vors_real,2);
            elseif  nargin ~= 1 || nargin ~= 5
                error('num of input arguments wrong');
            end
        end
        
        function [E k] = energySpectrum(obj)
            vorsAll = obj.vors_real+1i*obj.vors_imag;
            vorsCorrected = zeros(size(vorsAll));
            for i=1:obj.numFourier
                vorsCorrected(i,i:obj.numFourier,:,:) = vorsAll(i,1:obj.numFourier-i+1,:,:);
            end
            
            wavenumber = (1:size(vorsCorrected,2)-2)';
            waveFactor = wavenumber.*(wavenumber+1);
            E = zeros(size(wavenumber));
            
            for t=obj.t_start:obj.t_end
                for p=obj.p_lowb:obj.p_upperb
                    vors = reshape(vorsCorrected(2:size(vorsCorrected,1)-1,2:size(vorsCorrected,2)-1,p,t),...
                        size(vorsCorrected,1)-2,size(vorsCorrected,2)-2);
                    tmp = sum(abs(vors).^2,1)'./waveFactor;
                    E(1:length(E)) = E(1:length(E))+2*tmp(1:length(E));
                end
            end
            
            E = E/(obj.t_end-obj.t_start+1)/(obj.p_upperb-obj.p_lowb+1);
            k = wavenumber;
            obj.EperWavenumber = E;
            obj.wavenumber = wavenumber;
        end
        
        function nE = calEddyWavenumber(obj)
            waveNumF = obj.wavenumber.*(obj.wavenumber+1);
            n_n1 = sum(obj.EperWavenumber)./sum(obj.EperWavenumber./waveNumF);
            nE = 0.5*(sqrt(4*n_n1+1)-1);
            obj.energyContainWavenumber = nE;
        end
    end
end