站長留言

  • ✅ 本站維護及更新歷史紀錄,詳情請參考公告
  • ✅ 有任何意見、想法,歡迎留言給Spicy知道喔
  • ✅ 固定於每周一至周五更新Blogger文章,周末不定期
上課筆記高層大氣科學Matlab

【Matlab】GIM TEC 潮汐分解

tags: Matlab 高層大氣科學

Common 目標

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


GIM 資料 載點

潮汐分解方法 (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

沒有留言:

張貼留言

本網站建議使用電腦或平板瀏覽