【Matlab】GIM TEC 潮汐分解

tags: Matlab 高層大氣科學

Common 目標

1. 下載任一天GIM資料
   • 注意檔名結尾，應為 I.Z
   • 檔名範例：CODG0010.15I.Z
2. 進行潮汐分解 (UT)
3. 繪製 DW1 + SW2 + Xmean 每小時的變化圖，共25張)
4. 每張圖片上標示日期與時間

GIM 資料 載點

• 連結：ftp://ftp.aiub.unibe.ch/CODE

潮汐分解方法 (UT)

• DW1 (diurnal wave)：n=1, s=-1
• SW2 (semi-diurnal wave)：n=2, s=-2

副程式

data_read.m

%======================================
% Read the GIM data every one row line
%======================================
function [data_count,rdline]=data_read(data,fseg,data_count);
data_count=data_count+1;
line_index=(fseg(data_count)+1):(fseg(data_count+1)-1);
rdline='                                                                                                      ';
rdline(1:length(line_index))=data(line_index);

QRSol.m

function m_est=QRSol(A,B)

[Q,R]=qr(A);
sz=size(R);
U=R(1:sz(2),1:sz(2));
sz2=size(B);
for i=1:sz2(2)
   c=Q'*B(:,i);
   d=c(1:sz(2));
   m_est(:,i)=UTriSol(U,d);
end

  function x = UTriSol(U,b)
% x = UTriSol(U,b)
%
% Solves the nonsingular upper triangular system  Ux = b.
% where U is n-by-n, b is n-by-1, and X is n-by-1.

n = length(b);
x = zeros(n,1);
for j=n:-1:2
   x(j) = b(j)/U(j,j);
   b(1:j-1) = b(1:j-1) - x(j)*U(1:j-1,j);
end
x(1) = b(1)/U(1,1);

主程式，程式碼說明

1. 讀檔

%=================================================
% Read the GIM-TEC and do the tidal decomposition
%=================================================
clear all
close all

fn='CODG0010.15I';

fid=fopen(fn,'r');
data=fread(fid,inf);
fclose(fid);
fseg=find(data==10);
fseg=[0 fseg'];
data=char(data');
data_count=0;
headeryn=1;
count=1;

%----- Setting GIM parameters -----
gim_lat=87.5:-2.5:-87.5;
gim_lon=-180:5:180;

while data_count<(length(fseg)-1)
    [data_count,rdline]=data_read(data,fseg,data_count);
    %----- Read the GIM file header -----
    if headeryn==1
        if rdline(61:77)=='# OF MAPS IN FILE'
            map_num=str2num(rdline(1:6));
        end
        if rdline(61:73)=='END OF HEADER'
            headeryn=0;
        end
    else
        %----- Read GIM TEC data -----
        if rdline(61:71)=='END OF FILE'
            break
        end
        if rdline(61:76)=='START OF RMS MAP'
            break
        end
        if rdline(61:80)=='LAT/LON1/LON2/DLON/H'
            count2=1;
            for k=1:5
                [data_count,rdline]=data_read(data,fseg,data_count);
                tec_num=16;
                if k==5
                    tec_num=9;
                end
                for j=1:tec_num
                    index=1+(j-1)*5;
                    TEC(count,count2)=str2num(rdline(index:(index+4)))/10;
                    count2=count2+1;
                end
            end
            count=count+1;
        end
    end
end

2. 潮汐分解

%----- Setting the parameters of tidal decomposition -----
omega=2*pi/24;

%太陽同步與非太陽同步潮汐設定
n=1:1:3;
s=-4:1:4;
[n_matrix,s_matrix]=ndgrid(n,s);
param_len=length(n)*length(s);
n_matrix=reshape(n_matrix,[1 param_len]);
s_matrix=reshape(s_matrix,[1 param_len]);

%行星波(駐波)設定
n2=0;
s2=1:1:4;
[n_tmp,s_tmp]=ndgrid(n2,s2);
param_len2=length(n2)*length(s2);
n_tmp=reshape(n_tmp,[1 param_len2]);
s_tmp=reshape(s_tmp,[1 param_len2]);
n_matrix=[n_matrix n_tmp];
s_matrix=[s_matrix s_tmp];

param_len=param_len+param_len2;

%----- Calculate the GIM-TEC tidal composition -----
for i=1:length(gim_lat)
    lat=gim_lat(i);
    count=0;
    clear b
    Amatrix=zeros([length(gim_lon)*24 2*param_len+1]);
    for t=0:1:23
        index=i+t*length(gim_lat);
        for j=1:length(gim_lon)
            lon=gim_lon(j);
            count=count+1;
            b(count)=TEC(index,j);
            thi=n_matrix*omega*t-s_matrix*lon/180*pi;
            Amatrix(count,1:2:2*param_len)=cos(thi);
            Amatrix(count,2:2:2*param_len)=-1*sin(thi);
            Amatrix(count,2*param_len+1)=1;
        end
    end
    %左右兩邊乘上轉置矩陣，在左右對調，即可求得解
    x=QRSol(Amatrix'*Amatrix,Amatrix'*b');

3. 潮汐分解振幅與相位計算

    for k=1:param_len
        index=(k-1)*2+1;
        amp(i,k)=sqrt(x(index)^2+x(index+1)^2);
        theta(i,k)=atan(x(index+1)/x(index))/pi*180;
        if theta(i,k)>=0 & x(index+1)<0
            theta(i,k)=theta(i,k)-180;
        end
        if theta(i,k)<0 & x(index+1)>=0
            theta(i,k)=theta(i,k)+180;
        end
    end
    xmean(i)=x(end);
    
    clear Amatrix x b
end

4. 潮汐分解分量合併

• all tides 所有分量

% tide sum
for hour=0:24
    TEC2=zeros([length(gim_lat) length(gim_lon)]);  
    for ii=1:length(gim_lat)
        for jj=1:param_len
            TEC2(ii,:)=TEC2(ii,:)+amp(ii,jj)*cos(n_matrix(jj)*omega*hour-s_matrix(jj)*gim_lon*pi/180+theta(ii,jj)*pi/180);
        end
        TEC2(ii,:)=TEC2(ii,:)+xmean(ii);
    end

• 但實際上，目標要求：DW1 + SW2 + Xmean
  • 寫法1：用 find
    f=find(n_matrix == 1 & s_matrix == -1) | (n_matrix == 2 && s_matrix == -2)
  • 寫法2：用 if
    if (n_matrix(jj) == 1 && s_matrix(jj) == -1) || (n_matrix(jj) == 2 && s_matrix(jj) == -2)

% tide sum (DW1 + SW2 + xmean)
for hour=0:24
    TEC2=zeros([length(gim_lat) length(gim_lon)]);
    
    %DW1 (n=1, s=-1), SW2(n=2, s=-2)
    for ii=1:length(gim_lat)
        for jj=1:param_len
            if (n_matrix(jj) == 1 && s_matrix(jj) == -1) || (n_matrix(jj) == 2 && s_matrix(jj) == -2)
                TEC2(ii,:)=TEC2(ii,:)+amp(ii,jj)*cos(n_matrix(jj)*omega*hour-s_matrix(jj)*gim_lon*pi/180+theta(ii,jj)*pi/180);
            end
        end
        TEC2(ii,:)=TEC2(ii,:)+xmean(ii);
    end

5. 繪圖

    % fig
    [a,b]=contourf(gim_lon,gim_lat,TEC2,100);
    set(b,'edgecolor','none')
    caxis([0 60])
    colorbar
    colormap('jet')
    axis([-180 180 -90 90])
    set(gca,'xtick',-180:30:180,'ytick',-90:30:90,'fontsize',12)
    xlabel('Longitude (\circE)')
    ylabel('Latitude (\circN)')
    
    if hour<10
        hr=['0' num2str(hour)];
    else
        hr=num2str(hour);
    end
    
    title(['2015/01/01 ' hr '00UT'])
    %----- Plot the Global coast map -----
    load coast
    hold on
    plot(long,lat,'k')
    
    print('-dpng', ['GIM_' hr '.png'])
end

Result 結果

說明	動圖
DW1+ Xmean
SW2+ Xmean
DW1+ SW2+ Xmean (目標)
all tides
沒有 潮汐分解

Reference 參考資料

1. 高層大氣科學 上課講義
2. Atmospheric tide：https://en.wikipedia.org/wiki/Atmospheric_tide