Python Forum

xdata = [0.3063191028,-0.0156344344,-0.0155750443,-0.7206687321,-0.7473645659,-0.9174428618,-0.8839320182,-1.0399645639,-0.9997277955,-1.0157928079,-0.9888297188,-0.4533985964,0.0091163748,0.0026577054,0.5926386016,0.5992457462,1.004345373,0.9909529136,1.0392221881,1.0405287695,1.0606471537,1.0283835014,1.0149316519,0.9416591604,0.9685628594,0.9155869223,0.9088016075,0.6344640542,0.6142268898]
ydata = [1.1154790304,0.9978867036,1.0111900779,0.5702040049,0.5903372366,-0.0072010453,-0.0007720708,-0.45206232,-0.4390262009,-1.0375889614,-0.9978570085,-1.0612855969,-0.957932038,-0.904673998,-0.7489532501,-0.7689528542,-0.0266364437,-0.0265473586,0.2857701407,0.5784443763,0.5849624358,0.8579043226,0.8446900334,0.9316519346,0.9580805131,1.091470597,1.071560078,1.1199481327,1.0868233245]
xdata = np.array(xdata)
ydata = np.array(ydata)

def func1(coeff,x,y):
#     x = X[0]
#     y = X[1]

    n  = 8
#     % A  = ydata
#     % B  = -xdata
#     % C  = xdata. - ydata
#     % H  = zdata

    g = np.subtract(x,y)
    I_0 = np.subtract(x,y)   # x-y = C
    I_1 = np.multiply(I_0,coeff[2]) # c(x-y) = cC
    I_2 = np.multiply(coeff[1],-x)   #b(-x) = bB
    I_3 = np.multiply(coeff[0],y)  # aA


    I3_0 = np.subtract(I_1,I_2) # cC-bB
    I3_1 = np.subtract(I_3,I_1) # aA-cC
    I3_2 = np.subtract(I_2,I_3) # bB-aA

    I3_00 = np.multiply(I3_0,I3_1) # (cC-bB)(aA-cC)
    I3_01 = np.multiply(I3_00,I3_2) # (cC-bB)(aA-cC)(bB-aA)

    I3 = np.divide(I3_01,54) # (cC-bB)(aA-cC)(bB-aA)/54

    I2_0 = np.power((I3_1),2)  # (aA-cC)^2
    I2_1 = np.power((I3_0),2)  # (cC-bB)^2
    I2_2 = np.power((I3_2),2)  # (bB-aA)^2

    I2_00 = np.add(I2_0,I2_1)  # (aA-cC)^2 + (cC-bB)^2
    I2_01 = np.add(I2_00,I2_2) # (aA-cC)^2 + (cC-bB)^2 + (bB-aA)^2

    I2 = np.divide(I2_01,54)  # ((aA-cC)^2 + (cC-bB)^2 + (bB-aA)^2)/54

    th_0 = np.divide(I3,(np.power(I2,(3/2))))  # I3/(I2^(3/2))

#     print(th_0)

    th = np.arccos(np.clip((th_0),-1,1))  # arccos(I3/(I2^(3/2)))

#     print(th)

    ans_0 = np.divide(np.add((2*th),(np.pi)),6)   # (2*th + pi)/6
    ans_1 = np.divide(np.add((2*th),(3*np.pi)),6) # (2*th + 3*pi)/6
    ans_2 = np.divide(np.add((2*th),(5*np.pi)),6) # (2*th + 5*pi)/6

    ans_00 = np.multiply(np.cos(ans_0),2)  # 2*cos((2*th + pi)/6)
    ans_11 = np.multiply(np.cos(ans_1),2)  # 2*cos((2*th + 3*pi)/6)
    ans_22 = np.multiply(np.cos(ans_2),2)  # 2*cos((2*th + 5*pi)/6)

    ans_000 = np.power(np.absolute(ans_00),n)  # (abs(2*cos((2*th + pi)/6)))^n
    ans_111 = np.power(np.absolute(ans_11),n)  # (abs(2*cos((2*th + 3*pi)/6)))^n
    ans_222 = np.power(np.absolute(ans_22),n)  # (abs(2*cos((2*th + 5*pi)/6)))^n

    ans_0000 = np.add((np.power(np.absolute(ans_00),n)),(np.power(np.absolute(ans_11),n))) # (abs(2*cos((2*th + pi)/6)))^n + (abs(2*cos((2*th + 3*pi)/6)))^n 
    ans_1111 = np.add((ans_0000),(np.power(np.absolute(ans_22),n)))  # (abs(2*cos((2*th + pi)/6)))^n + (abs(2*cos((2*th + 3*pi)/6)))^n + (abs(2*cos((2*th + 5*pi)/6)))^n

    sna_0 = np.power(np.multiply(3,I2),(n/2))  # (3*I2)^(n/2) !!
    sna_1 = 2*(np.power(1.0167,n)) # 2*(sigma^n) !!

    sna_00 = np.multiply(sna_0,ans_1111)
    sna_11 = np.subtract(sna_00,sna_1)

    return sna_11
x0 = np.array([1.0, 1.0, 1.0])
res_lsq = least_squares(func1, x0,loss='cauchy',f_scale=0.001,args=(xdata, ydata))
res_lsq.x

xmin, xmax, ymin, ymax = xdata.min(), xdata.max(), ydata.min(), ydata.max()
g = lambda x, y: func1(res_lsq.x, x, y)**2
X, Y = np.meshgrid(np.linspace(xmin-0.2, xmax+0.2, 1000), np.linspace(ymin-0.2, ymax+0.2, 1000))
Z = g(X, Y)  # heavy computations.... 

# plotting the result 
plt.contourf(X, Y, np.log(Z), levels=[-10,-8])
plt.gca().scatter(xdata, ydata)
plt.colorbar()
plt.show()