import random
import math
import copy
import time
import numpy as np
import scipy
import matplotlib.pyplot as plt

global profile
profile =[]
global dx 
dx =0 

def func( i, x, f, df ):
    temp1 = -(1/2/math.pi**2) * math.sin(f[i])**2 * df[i]  * 4* math.pi
    temp2 = -(1/2/math.pi**2) * math.sin(f[i+1])**2 * df[i+1] * 4* math.pi
    # print("temp1 = %f, temp2 = %f "%(temp1,temp2))
    temp = (temp1+temp2)*(x[i+1] - x[i])/2
    return temp


def baryonNumber(x, f):
    i=0
    temp = 0
    for xi in x:
        if i==len(x)-1:
            break
        temp = temp + f(i, x, profile, df)
        i=i+1
    return temp

r0 = 1e-5
r1 = 10
N = 100
x = np.linspace(r0, r1, N)
profile = np.zeros_like(x)
df = np.zeros_like(x)
i=0
for xi in x:
    profile[i] = math.pi/(r1-r0) * ( r1 - x[i])
    df[i] = -math.pi/(r1-r0)
    i=i+1


baryon = baryonNumber(x, func)
print(baryon)
i=5
while True:
    pi = profile[i]
    profile[i] = profile[i] + random.uniform(-1,1) * 0.5
    temp = baryonNumber(x, func)
    if abs(temp-baryon)<1e-3:
        print("原来为 %f, 更新后的解为 %f "%(pi, profile[i]) )
        break
    else:
        print("结果不符, 现在值为 %f"%( profile[i] ) )
        profile[i] = pi



print(baryon)