关于用Fujita-Takahashi公式建立台风风场模型的问题

GGGrADS

GGGrADS 发表于 2020-07-29 16:05
纬度转换为弧度是用lat*pi/180吗

我在matlab没找到np.sin

来自安卓APP客户端

GGGrADS

islandowner 发表于 2020-07-10 11:33
f = 2.0*omega*np.sin(Lat)
omega     = 7.292E-05   
Lat是弧度不是角度

距离r和最大风速半径R用km,
Vx,Vy用m/s,
x和y,x0和y0我用的utm坐标转换成了m，
气压用hPa
空气密度用的kg/m3
这些量的单位是这样吗，我算不对。。。

来自安卓APP客户端

islandowner

GGGrADS 发表于 2020-7-29 16:05
纬度转换为弧度是用lat*pi/180吗

是的，角度转弧度是这样，不过你要确定matlab的sin函数输入的是弧度还是角度

islandowner

GGGrADS 发表于 2020-7-29 16:08
我在matlab没找到np.sin

np.sin是python语言的，其实就是sin三角函数

islandowner

GGGrADS 发表于 2020-7-29 16:11
距离r和最大风速半径R用km,
Vx,Vy用m/s,
x和y,x0和y0我用的utm坐标转换成了m，

r和R用m和km无所谓，因为他们相除无量纲了。Vx Vy用m。气压用Pa，密度用kg/m3。然后公式里面乘10的三次幂不要。

GGGrADS

islandowner 发表于 2020-07-29 22:18
r和R用m和km无所谓，因为他们相除无量纲了。Vx Vy用m。气压用Pa，密度用kg/m3。然后公式里面乘10的三次幂不要。

但是在Ω那一项不是有r*R和R平方的项吗，气压用pa不用hPa吗，为啥Ω公式里10的三次幂要去掉

来自安卓APP客户端

GGGrADS

islandowner 发表于 2020-7-29 22:18
r和R用m和km无所谓，因为他们相除无量纲了。Vx Vy用m。气压用Pa，密度用kg/m3。然后公式里面乘10的三次幂 ...

clear all
CMATC=read_CMABST('E:\X_file\tf\cntf',2018);
lon=CMATC{21,1}(:,8);
lat=CMATC{21,1}(:,7);
wind=CMATC{21,1}(:,10);
hPa=CMATC{21,1}(:,9);
ommiga=7.292e-5;

varargin=[lat lon];
[x,y,z]=ll2utm(varargin);
x1=x;
y1=y;
[oo,pp]=size(varargin);
r1=matrix_hav(varargin)*1000;
r2=[r1 r1(end)];
r=reshape(r2,[length(r2),1]);
v=r/(6*60*60);

seigm1=matrix_fangweijiao(varargin);
seigm2=[seigm1 seigm1(end)];
seigm=reshape(seigm2,[length(seigm2),1]);

A=[v seigm];
[vx,vy]=V_XY(A);

ra_lat=radian2degree(lat);

latmin=20;
latmax=33;
lonmin=117;
lonmax=133;
step=1;
s2=1/step;
LAT=latmin.*s2-1;
LON=lonmin.*s2-1;
x_lat=(latmin:step:latmax);
y_lon=(lonmin:step:lonmax);
C1=1;
C2=0.8;
density=1.293;
theta=20;

% R=28.52*tanh(0.0873*(lon-28))+12.22*exp((hPa-1013.2)/33.86)+38.2;
R=1119*(1010-hPa).^(-0.805);

for j=1:length(lon);
   for r=latmin:step:latmax;
      for t=lonmin:step:lonmax;
        distance(ceil(r.*s2-LAT),ceil(t.*s2-LON),j)=hav(lat(j),lon(j),r,t);
      end
   end
end


for m=1:length(lon);
    m;
    distance_a=distance(:,:,m);%  r
    [srow scol]=find(distance_a<=2*R(m));%s=short,l=long
    [lrow lcol]=find(distance_a>2*R(m));
    dw=wind(m);
    dR=R(m);%  R
    Vx=vx(m);
    Vy=vy(m);
    x0=x1(m);
    y0=y1(m);
    P=1013.25-hPa(m);
   
   
    for i=1:length(srow);
        i;
        xx=srow(i);
        yy=scol(i);
        var=[x_lat(xx) y_lon(yy)];
        Ff=2*ommiga*sin(radian2degree(x_lat(xx)));
        [x2 y2 z]=ll2utm(var);
        r=distance_a(xx,yy);
        Rr=abs(r-dR);
        omiga1=-0.5*Ff+sqrt(0.25*Ff^2+1000*2*P*(1+2*(r^2/dR^2))^(-1.5)/(density*dR^2));
        Wx(xx,yy,m)=C1*Vx*exp(-0.25*pi*Rr/dR)-C2*omiga1*((x2-x0)*sind(theta)+(y2-y0)*cosd(theta));
        Wy(xx,yy,m)=C1*Vy*exp(-0.25*pi*Rr/dR)-C2*omiga1*((x2-x0)*cosd(theta)-(y2-y0)*sind(theta));
%         W(ceil(r.*s2-LAT),ceil(t.*s2-LON))=sqrt(Wx(r.*s2-LAT,t.*s2-LON)^2+Wy(r.*s2-LAT,t.*s2-LON)^2);
    end
    for i=1:length(lrow);
        i;
        xx=lrow(i);
        yy=lcol(i);
        var=[x_lat(xx) y_lon(yy)];
        Ff=2*ommiga*radian2degree(x_lat(xx));
        [x2 y2 z]=ll2utm(var);
        r=distance_a(xx,yy);
        Rr=abs(r-dR);
        omiga2=-0.5*Ff+sqrt(0.25*Ff^2+1000*P/(density*r*dR*(1+r/dR)^2));
        Wx(xx,yy,m)=C1*Vx*exp(-0.25*pi*Rr/dR)-C2*omiga2*((x2-x0)*sind(theta)+(y2-y0)*cosd(theta));
        Wy(xx,yy,m)=C1*Vy*exp(-0.25*pi*Rr/dR)-C2*omiga2*((x2-x0)*cosd(theta)-(y2-y0)*sind(theta));
    end
end

我这样算出来风速上千。。。