风场反演python代码

import math
import numpy as np
import operator
from itertools import chain
import matplotlib.pyplot as plt
def degrees_to_radians(degrees):
    radians = degrees * math.pi / 180
    return radians




if __name__ == '__main__':
    #x = np.array( [[degrees_to_radians(0)], [degrees_to_radians(90)], [degrees_to_radians(180)], [degrees_to_radians(270)], [degrees_to_radians(360)]])

    #x1 = np.array([[degrees_to_radians(135)], [degrees_to_radians(90)], [degrees_to_radians(0)], [degrees_to_radians(45)],[degrees_to_radians(180)]])
    # x = np.array([[degrees_to_radians(0)], [degrees_to_radians(45)], [degrees_to_radians(90)], [degrees_to_radians(135)],
    #      [degrees_to_radians(180)]])
    # x = np.array([[math.radians(0)], [math.radians(45)], [math.radians(90)], [math.radians(135)],   #math.radians()将角度转换为弧度
    #      [math.radians(180)]])
    #print(x[1])
    #y = np.mat([[-0.188288, -0.176802], [0.825875, 0.857461], [-0.248817, -0.521539], [-1.12957, -0.989984],[-0.176802,0.171071]])
    #y = np.mat([[-0.594857], [0.490327], [-0.484808], [-1.73742]])  #ID1测试，vh=3.26,vf=-0.61,vd=β0=92.82
    #y1 = np.mat([[0.86], [1.78], [0.21], [0.85],[-1.01]])
    #y1 = np.mat([[0.21], [0.85], [1.78], [0.12], [-1.01]])
    #x = np.array([[degrees_to_radians(0)], [degrees_to_radians(90)], [degrees_to_radians(180)], [degrees_to_radians(270)]])
    # #y = np.mat([[-1.27], [-1.06], [-1.44],[-2.53]]) #第一列
    # # #y = np.mat([[-1.95], [0.86], [0.99], [-2.6]])
    # y = np.mat([[-1.01], [-6.42], [-2.14], [3.21]])  #最后一列


    #y = np.mat([[-2.07], [-3.63], [-0.84], [1.53]])

    x = np.array(
          [[degrees_to_radians(0)], [degrees_to_radians(45)],[degrees_to_radians(90)], [degrees_to_radians(135)],[degrees_to_radians(180)], [degrees_to_radians(225)],[degrees_to_radians(270)],[degrees_to_radians(315)]])
    # x = np.array(
    #     [[degrees_to_radians(0)], [degrees_to_radians(44.8)],[degrees_to_radians(89.9)], [degrees_to_radians(134.9)],[degrees_to_radians(179.9)], [degrees_to_radians(224.8)],[degrees_to_radians(269.9)],[degrees_to_radians(314.9)]])
    # # y = np.mat([[-0.4], [0.37], [0.53], [-0.33],[-1.47],[-2.21],[-2.6],[-1.6]])
    y = np.mat(
        [[0.122175956 + 10], [-0.396681876 + 10], [-0.726924423 + 10], [-0.917665079 + 10], [0.088133387 + 10],
         [0.745351522 + 10], [0.872576389 + 10], [0.669659167 +10]])
    # y = np.mat([[-2.1073], [3.0352], [2.6076], [-2.1743]])
    y1 = np.mat([[0.122175956+0.7], [-0.396681876+0.7], [-0.726924423+0.7], [-0.917665079+0.7], [0.088133387+0.7], [0.745351522+0.7], [0.872576389+0.7], [0.669659167+0.7]])
    # y = np.mat([[-2.1073], [3.0352], [2.6076], [-2.1743]])
    # y = np.mat([[0.604227], [0.579773], [-0.59057], [-0.51622]])
    #y = np.mat([[0.195], [0.685], [-0.823], [-0.15]])
    A = np.mat([math.cos(x[0]), math.sin(x[0]), 1])
    B = 0
    for i in range(1, len(x)):
        B = np.mat([math.cos(x[i]), math.sin(x[i]), 1])
        A = np.vstack((A, B))  # 矩阵纵向拼接
    print(A)
    print("-------------------------------------")
    L = y[:, 0]    # 第1列数据
    L1 = y1[:, 0]
    #A_in= np.linalg.inv(A)
    W = np.dot((np.linalg.inv(A.T * A)), (A.T))
    W1 = W*L
    W2 = W*L1
    p = np.dot((np.linalg.inv(A.T * A)), (A.T * L))
    #p = np.dot((np.linalg.inv(A.T * A)), (A.T * L))
    #y = np.round(y, decimals=8)
    print(p)
    a1 = math.sqrt(p[0]*p[0]+p[1]*p[1])
    # # 使用 math.atan2 计算弧度值（范围从 -π 到 π）,math.atan2函数计算角度，因为它可以处理所有四个象限的情况，并返回正确的角度
    # xc = math.atan2(-p[0] ,p[1])    # xc弧度
    # # 将弧度转换为度（范围从 -180 到 180）
    # xc = math.degrees(xc)
    # # 如果结果小于0（即第二或第三象限），则加上360度以得到0-360度的范围
    # if xc < 0:
    #     xc += 360

        #反正切函数将xc限制在了（-Π/2，Π/2）内，但tan(xc+n*Π)=tan(xc)=-p[0] / p[1]
    if ((-p[0] / p[1])>0 and -p[0]>0): #(p[1],-p[0])在第一象限
        xc = math.atan(-p[0] / p[1])
    elif ((-p[0] / p[1])>0 and -p[0]<0):  #math.atan(-p[0] / p[1])>0,(p[1],-p[0])在第三象限，
        xc = math.atan(-p[0] / p[1]) + math.pi
    elif ((-p[0] / p[1])<0 and -p[0]>0):  #math.atan(-p[0] / p[1])<0,(p[1],-p[0])在第二象限
        xc = math.atan(-p[0] / p[1])+math.pi
    elif ((-p[0] / p[1])<0 and -p[0]<0):  #(p[1],-p[0])在第四象限
        xc = math.atan(-p[0] / p[1])

    if (xc  + math.pi / 2) > 0 and (xc  + math.pi / 2) < 2 * math.pi:
        β0  = xc  + math.pi/ 2
    elif (xc  + math.pi/ 2) > 2 * math.pi:
        β0  = xc  + math.pi/ 2 - 2 * math.pi;
    elif ((xc  + math.pi / 2) < 0):
        β0  = xc  + math.pi/ 2 + 2 * math.pi;  #β0弧度

    y0 = float(p[2])
    β0 = float(β0)
    #print("y0=",y0)
    vf = y0/math.cos(degrees_to_radians(20))   # 仰角20，垂直风速度
    vh = a1/math.sin(degrees_to_radians(20))   # 水平风速
    vd = 180*β0/math.pi   # 水平风向β0转角度
    #vd = math.degrees(b2)  # 水平风向β0角度

    print("p[0]:")
    print(p[0])
    print("p[1]:")
    print(p[1])
    print("xc:")
    print(math.degrees(xc))
    print("振幅A，水平风速:")
    print(a1)
    print(vh)
    print("相位b，即水平风向:")
    print(vd)
    print("截距y0,垂直风速:")
    print(y0)
    print(vf)

    # 绘制原始数据点
    # x2 = [x[0], x[1], x[2], x[3]]
    # y2 = [y[0][0], y[1][0], y[2][0], y[3][0]]
    x2 = [x[0], x[1], x[2], x[3], x[4], x[5], x[6], x[7]]
    y2 = [y[0][0], y[1][0], y[2][0], y[3][0],y[4][0], y[5][0], y[6][0], y[7][0]]

    plt.scatter(x2, y2, color='red', label='Original Data')

    # 绘制拟合的余弦曲线
    #x_fit = np.array(
      #  [[degrees_to_radians(0)], [degrees_to_radians(90)], [degrees_to_radians(180)], [degrees_to_radians(270)]])
    x_fit = np.linspace(0, 2 * np.pi, 50)

    y_fit = y0+a1 * np.sin(x_fit - xc)

    plt.plot(x_fit, y_fit, color='blue', label='Fitted Curve')

    plt.xlabel('Degrees')
    plt.ylabel('Values')
    plt.legend()
    plt.show()