- 帖子
- 4
- 精华
- 0
- 积分
- 40
- 阅读权限
- 10
- 注册时间
- 2017-6-22
- 最后登录
- 2017-8-23
|
本帖最后由 方自之 于 2017-6-28 19:00 编辑
各位老铁、Crossin先生好:
py新手,由于之前都是用fortran做数值模拟较多,有些习惯延续到python.我的问题如下:
1、如代码,我遍历数组p[3000,2000]依然是fortran的习惯是:两个嵌套循环。python中是否建议这样使用,是否有更规范便捷的方法?
2、如代码,我在两个for循环内的if语句,报语法错误,很捉急,没有看出来,请大家帮忙找茬,第35行。
3、由于数据较大,我起初想创建的数组维度是p[240,3000,2000],但报错memory error,自己查了(其实根据字面意思也能猜出来)是数组太大,我无奈改为p[3000,2000]。请问python中该如何解决?
4、代码很简单,都是一元函数运算,如果大家觉得在书写方法上有可以优化的地方,或者有些坏习惯,也请告诉我。
先谢谢各位了- # -*- coding: utf-8 -*-
- #本程序用于实现理想台风模拟,参考防灾减灾手册。
- import numpy as np
- import matplotlib.pyplot as plt
- import matplotlib.mlab as mlab
- from mpl_toolkits.basemap import Basemap, cm
- import math
- #参数设定
- R_max = 39.0 #最大风速半径(km)
- pi = 3.1415926535897932384626 #pai
- P_max = 1010.0 #无穷远处最大气压(hPa)
- P_0 = 965.0 #台风中心气压
- V_0 = 25.0 #台风移速(km/h)
- x_def = 2000
- y_def = 3000 #区域范围(km)
- t_def = 120 #模拟时间区间(h)
- theta_1 = pi/18.0
- theta_2 = pi/7.2 #流入角度
- W_r = 3.029*(P_max-P_0)**0.644
- #Surface Pressure(hPa)
- P = np.zeros([y_def,x_def])
- U = np.zeros([y_def,x_def])
- V = np.zeros([y_def,x_def])
- pit_x = 1000
- pit_y = 1500
- P[pit_y,pit_x] = P_0 #设置台风中心
- for i in range(x_def):
- for j in range(y_def):
- r=math.sqrt(abs(i-pit_x)**2+abs(j-pit_y)**2)
- theta_3 = pi/12.0*((r-R_max)/(0.2*R_max)+(pi/18)
- if r>0 and r<= R_max:
- A = -((i-pit_x)*math.sin(theta_1)+(j-pit_y)*math.cos(theta_1))
- B =(i-pit_x)*math.cos(theta_1)-(j-pit_y)*math.sin(theta_1)
- P[j,i] = P_0+(0.25*(P_max-P_0)*(r/R_max)**3)
- U[j,i] = (R_max/(r-R_max)0+W_r*(r/R_max)**1.5*A/r
- V[j,i] = (R_max/(r-R_max)*V_0)+W_r*(r/R_max)**0.5*B/r
- elif(r>R_max and r<R_max*1.2): #theta_3
- P[j,i] =P_max-(0.75*(P_max-P_0)*(R_max/r))
- A = -((i-pit_x)*math.sin(theta_3)+(j-pit_y)*math.cos(theta_3))
- B =(i-pit_x)*math.cos(theta_3)-(j-pit_y)*math.sin(theta_3)
- U[j,i] = (R_max/(r-R_max)0+W_r*(r/R_max)**1.5*A/r
- V[j,i] = (R_max/(r-R_max)*V_0)+W_r*(r/R_max)**0.5*B/r
- elif(r>=1.2*R_max): #theta_2
- P[j,i] =P_max-(0.75*(P_max-P_0)*(R_max/r))
- A = -((i-pit_x)*math.sin(theta_2)+(j-pit_y)*math.cos(theta_2))
- B =(i-pit_x)*math.cos(theta_2)-(j-pit_y)*math.sin(theta_2)
- U[j,i] = (R_max/(r-R_max)0+W_r*(r/R_max)**1.5*A/r
- V[j,i] = (R_max/(r-R_max)*V_0)+W_r*(r/R_max)**0.5*B/r
- #Print Projection
复制代码 |
|