#! /usr/bin/env python3
# -*- coding: utf-8 -*-
# vim:fenc=utf-8
#
# Author: 冯振华
# StudentNumber: 202221140015
# Verdion: V4.0
# Copyright © 2023 feng <feng@arch>
#
# Distributed under terms of the MIT license.

import numpy as np
import matplotlib.pyplot as plt

# Define the parameters
N = 20 # number of spins in each direction
J = -1 # coupling constant
beta = 1 # inverse temperature
nsteps = 1000 # number of Monte Carlo steps

# Initialize the spins randomly
spins = np.random.choice([-1, 1], size=(N, N))

# Define the energy function
def energy(spins):
    return -J * (np.sum(spins[:-1,:] * spins[1:,:]) + np.sum(spins[:,:-1] * spins[:,1:])) \
           - J * (np.sum(spins[0,:] * spins[-1,:]) + np.sum(spins[:,0] * spins[:,-1]))

# Define the Metropolis algorithm
def metropolis(spins, beta):
    for i in range(N):
        for j in range(N):
            # Choose a random spin to flip
            x, y = np.random.randint(0, N), np.random.randint(0, N)
            # Calculate the energy difference
            delta_E = 2 * J * spins[x, y] * \
                      (spins[(x+1)%N, y] + spins[(x-1)%N, y] + spins[x, (y+1)%N] + spins[x, (y-1)%N])
            # Flip the spin with probability according to the Boltzmann factor
            if np.random.uniform(0, 1) < np.exp(-beta * delta_E):
                spins[x, y] = -spins[x, y]
    return spins

# Run the Monte Carlo simulation
energies = []
for i in range(nsteps):
    spins = metropolis(spins, beta)
    energies.append(energy(spins))

# Calculate the average energy
E = np.mean(energies)

# Print the result
print(f"Average energy: {E:.4f}")

# Plot the energy as a function of Monte Carlo steps
plt.plot(energies)
plt.xlabel("Monte Carlo step")
plt.ylabel("Energy")
plt.show()

# Plot the spin configuration
plt.imshow(spins, cmap='binary')
plt.show()

#! /usr/bin/env python3
# -*- coding: utf-8 -*-
# vim:fenc=utf-8
#
# Copyright © 2023 feng <feng@archlinux>
#
# Distributed under terms of the MIT license.
"""
作者：冯振华
学号：202221140015
版本：V4.1
日期：2023年05月20日
项目：Monte Carlo methods for Overcoming critical slowing down
"""
import numpy as np
import matplotlib.pyplot as plt

def initialize_spins(N):
    """
    Initializes the spin lattice with random spins of -1 or +1
    """
    return np.random.choice([-1, 1], size=(N, N))

def bond_probability(beta):
    """
    Calculates probability of forming a bond between two neighboring spins
    """
    return 1 - np.exp(-2*beta)

def update_spins(spins, beta):
    """
    a single update of the spins using the Wolff algorithm
    """
    N = spins.shape[0]
    visited = np.zeros((N, N), dtype=bool)
    seed = (np.random.randint(N), np.random.randint(N))
    cluster_spins = [seed]
    cluster_perimeter = [(seed[0], (seed[1]+1)%N), (seed[0], (seed[1]-1)%N),
                         ((seed[0]+1)%N, seed[1]), ((seed[0]-1)%N, seed[1])]
    while len(cluster_perimeter) > 0:
        i, j = cluster_perimeter.pop()
        if not visited[i,j] and np.random.random() < bond_probability(beta):
            cluster_spins.append((i, j))
            visited[i, j] = True
            cluster_perimeter += [(i, (j+1)%N), (i, (j-1)%N),
                                  ((i+1)%N, j), ((i-1)%N, j)]
    spins[tuple(zip(*cluster_spins))] *= -1
    return spins

def calculate_energy(spins,J):
    """
    Calculates the energy of the spin configuration
    """
    N = spins.shape[0]
    energy = 0
    for i in range(N):
        for j in range(N):
            spin = spins[i, j]
            neighbors = spins[(i+1)%N, j] + spins[i, (j+1)%N] + spins[(i-1)%N, j] + spins[i, (j-1)%N]
            energy += 2*J*spin * neighbors    
    return energy 

def run_simulation(N=20, J=-1,beta=1.0, num_steps=2000, num_measurements=100):
    """
    Performs a full simulation of the Wolff algorithm on the spin lattice
    """
    spins = initialize_spins(N)
    energies = []
    mags = []
    for step in range(num_steps):
        spins = update_spins(spins, beta)
        if step % (num_steps // num_measurements) == 0:
            energy = calculate_energy(spins,J)
            mag = np.sum(spins)
            energies.append(energy)
            mags.append(mag)
    return energies, mags 

energies, mags = run_simulation(N=20,J=-1,beta=1.0, num_steps=100000, num_measurements=100)
#
# Calculate the average energy
E = np.mean(energies)
# Print the result
print(f"Average energy: {E:.4f}")
# Plot the energy as a function of Monte Carlo steps
plt.plot(energies)
plt.xlabel("Monte Carlo step")
plt.ylabel("Energy")
plt.show()

# Plot the spin configuration
plt.plot(mags)
plt.xlabel("Monte Carlo step")
plt.ylabel("Spin")
plt.show()
"""
    To run the simulation, simply the `run_simulation` function with the desired parameters, such as `energies, mags = run_simulation(N=30,=1.0, num_steps=100000, num_measurements=100)`. This will simulate 100,000 Wolff updates on a 30x30 lattice at a temperature of 1. and take measurements of the energy and magnetization every 1,000 steps. The output is two arrays, `energies` and `mags`, that record the energy and magnetization at each measurement step. These arrays can be used to averages and variances to analyze the system properties.
"""

#! /usr/bin/env python3
# -*- coding: utf-8 -*-
# vim:fenc=utf-8
# Author:        冯振华 & 程国栋 
# StudentNumber: 202221140015 & 202221140011
# Version:       v2.2  for linux
# Linux:         冯振华
# Windows:       程国栋 
# Date:          2023年05月18日星期四多云北京市北京师范大学
# Copyright © 2023 feng <feng@archlinux>
#
# !!注意：本脚本在linux编写，vim编辑器下
#       :set ff      #查看文件格式
#       :set ff=dos  #设置为dos格式
#       :set ff=unix #设置为unix格式
# 如果不设置正确的格式地，运行程序会报错
#
import math
import matplotlib.pyplot as plt
import os
import requests
import json
# 输入初速度大小和空气阻力系数
v_0 = input("请输入初速度大小: ")
v_0 = float(v_0)
k = input("请输入空气阻力系数: ")
k = float(k)
# 获取本机Ip对应的地理位置
Ip_data=os.popen('curl -s  cip.cc').readlines()
Ip_local=Ip_data[1].split(' ')
city=Ip_local[-1].strip()
# 从高德地图获取纬度
GaoDeKey = '79e3576b80e0b544d0d17dd643e4654e'
url = 'https://restapi.amap.com/v3/geocode/geo'
params = {'key': GaoDeKey,
    'address': city}
res = requests.get(url,params)
jd = json.loads(res.text)
coords = jd['geocodes'][0]['location']
jd_value = coords.split(',')
latitude=float(jd_value[1])
latitude = latitude/180*math.pi
# 根据纬度计算本地的重力加速度g
g=9.7803*(1+0.005324*((math.sin(latitude))**2)-0.000005*((math.sin(2*latitude))**2))
# 计算最大水平位移和初速度与水平方向夹角的关系
def x_max(theta):
    rad = math.radians(theta)
    return (v_0**2* math.sin(2*rad)/ g)
    angles = range(0, 91, 1)
    x = [x_max(angle) for angle in angles]
#with the drag force
def simulate(v_1, theta,y_height):
    theta = math.radians(theta)
    v_x = v_1 * math.cos(theta)
    v_y = v_1 * math.sin(theta)
    x = 0
    y = float(y_height)
    t = 0
    delta_t = 0.01
    traj = [(x, y)]
    while y >= 0:
        a_x = -k * v_x ** 2 / (2 * v_1)
        a_y = -g - k * v_y ** 2 / (2 * v_1)
        v_x += a_x * delta_t
        v_y += a_y * delta_t
        x += v_x * delta_t
        y += v_y * delta_t
        t += delta_t
        traj.append((x, y))
    return traj
# 将各种情况集中于一张图上
plt.figure(figsize=(200,200),dpi=100)
plt.subplots_adjust(left=None,bottom=None,right=None,top=None,wspace=None,hspace=0.5)
plt.figure(1)
ax1 = plt.subplot(221)
ax1.margins(0.15)
ax1.plot(angles, x)
ax1.set_title("$\\theta$ with $x_{max}$ Free")
ax1.set_xlabel("$\\theta$(degrees)")
ax1.set_ylabel("$x_{max}$(meters)")
ax2 = plt.subplot(222)
ax2.margins(0.15)
angles = range(0, 91, 10)
for theta in angles:
    traj = simulate(v_0, theta,0) # 地表
    x_s = [t[0] for t in traj]
    y_s = [t[1] for t in traj]
    ax2.plot(x_s, y_s, label=f"${theta}^\circ$")
ax2.legend(loc="upper left")
ax2.set_title(f"Projectile Motion with Air Resistance ($v_0={v_0} m/s$)")
ax2.set_xlabel("Distance (meters)")
ax2.set_ylabel("Height (meters)")
ax3 = plt.subplot(223)  
ax3.margins(0.15)
angles = range(0, 91, 10)
for theta in angles:
    traj = simulate(v_0, theta,1)  # 100m高塔
    x_s = [t[0] for t in traj]
    y_s = [t[1] for t in traj]
    ax3.plot(x_s, y_s, label=f"${theta}^\circ$")
ax3.legend(loc="upper left")
ax3.set_title(f"Projectile Motion with Air Resistance ($v_0={v_0} m/s$)")
ax3.set_xlabel("Distance (meters)")
ax3.set_ylabel("Height (meters)")
ax4 = plt.subplot(224)
ax4.margins(0.15)
angles = range(0, 91, 10)
for theta in angles:
    traj = simulate(v_0, theta,2)  # 200m高塔
    x_s = [t[0] for t in traj]
    y_s = [t[1] for t in traj]
    ax4.plot(x_s, y_s, label=f"${theta}^\circ$")
ax4.legend(loc="upper left")
ax4.set_title(f"Projectile Motion with Air Resistance ($v_0={v_0} m/s$)")
ax4.set_xlabel("Distance (meters)")
ax4.set_ylabel("Height (meters)")
plt.show()

#! /usr/bin/env python3
# -*- coding: utf-8 -*-
# vim:fenc=utf-8
# Author:        冯振华 & 程国栋 
# StudentNumber: 202221140015 & 202221140011
# Version:       v2.2  for windows
# Linux:         冯振华
# Windows:       程国栋 
# Date:          2023年05月18日星期四多云北京市北京师范大学
# Copyright © 2023 feng <feng@archlinux>
#
# !!注意：本脚本在linux编写，vim编辑器下
#       :set ff      #查看文件格式
#       :set ff=dos  #设置为dos格式
#       :set ff=unix #设置为unix格式
# 如果不设置正确的格式地，运行程序会报错
#
import math
import numpy as np
import matplotlib.pyplot as plt
import requests
import json
from bs4 import BeautifulSoup
from scipy.integrate import odeint


# 输入初速度大小和空气阻力系数
v_0 = float(input("请输入初速度大小: "))
k = float(input("请输入空气阻力系数: "))


# 获取本机Ip对应的地理位置
url = 'http://www.cip.cc/'
header = {'User-Agent': 'Mozilla/5.0 (Windows NT 10.0; Win64; x64; rv:109.0) Gecko/20100101 Firefox/112.0'}
html = requests.get(url, headers=header)  # 访问网址
bsObj = BeautifulSoup(html.content, 'lxml')  # 使用lxml解析器
msg = bsObj.find('div', class_="data kq-well")
msg = [i for i in msg]
city = msg[1]


# 从高德地图获取纬度
GaoDeKey = '79e3576b80e0b544d0d17dd643e4654e'
url = 'https://restapi.amap.com/v3/geocode/geo'
params = {'key': GaoDeKey,
          'address': city}
res = requests.get(url, params)
jd = json.loads(res.text)
coords = jd['geocodes'][0]['location']
jd_value = coords.split(',')
latitude = float(jd_value[1])
latitude = latitude/180*math.pi  # 转换成弧度


# 根据纬度计算本地的重力加速度g
g = 9.7803*(1+0.005324*((math.sin(latitude))**2)-0.000005*((math.sin(2*latitude))**2))


def simulate(v_1, theta, y_height):
    # 轨迹图
    theta = math.radians(theta)
    v_x0 = v_1 * math.cos(theta)
    v_y0 = v_1 * math.sin(theta)

    def motion1(w, t, k):
        v_x, v_y = w
        return[-k * v_x ** 2, -g - k * v_y ** 2]
    t = np.linspace(0, 20, 300000)
    track1 = odeint(motion1, (v_x0, v_y0), t, args=(k, ))
    v_x1 = [i for i in track1[:, 0]]
    v_y1 = [j for j in track1[:, 1] if j > 0]
    v_x1 = v_x1[:len(v_y1)]
    delta_t = (20 - 0) / 300000
    i = 1
    x = 0
    y = y_height
    y_list1 = [y]
    x_list1 = [0]
    while i < len(v_y1) - 1:
        x += ((v_x1[i] + v_x1[i-1]) / 2) * delta_t
        y += ((v_y1[i] + v_y1[i-1]) / 2) * delta_t
        if y < 0:
            break
        x_list1.append(x)
        y_list1.append(y)
        i += 1

    def motion2(w, t, k):
        v_x, v_y = w
        return[-k * v_x ** 2, -g + k * v_y ** 2]

    track2 = odeint(motion2, (v_x1[-1], 0), t, args=(k, ))
    v_x2 = [i for i in track2[:, 0]]
    v_y2 = [j for j in track2[:, 1] if j < 0]
    j = 1
    x = x_list1[-1]
    y = y_list1[-1]
    y_list2 = [y]
    x_list2 = [x]
    while j < len(v_y2) - 1:
        x += ((v_x2[j] + v_x2[j-1]) / 2) * delta_t
        y += ((v_y2[j] + v_y2[j-1]) / 2) * delta_t
        if y < 0:
            break
        x_list2.append(x)
        y_list2.append(y)
        j += 1
    x_list = x_list1 + x_list2
    y_list = y_list1 + y_list2
    return x_list, y_list


def x_MaxTheta(v, height):
    x_max = []
    for angle in range(1, 91, 1):
        x_list, y_list = simulate(v_1=v, theta=angle, y_height=height)
        x_max.append(x_list[-1])
    return x_max


heght = float(input('输入初始高度:'))
x_max = x_MaxTheta(v=v_0, height=heght)
plt.plot(range(1, 91, 1), x_max, 'g-')
plt.show()
x = 0
j = 1
for i in x_max:
    j += 1
    if i > x:
        x = i
        num = j - 1
print('当前初始条件下抛得最远的角度是{},最远距离是{}'.format(num, x))
print('输入要模拟的轨迹的抛出角度:')
a = float(input())
x, y = simulate(v_1=v_0, theta=a, y_height=heght)
plt.plot(x, y, 'r-')
plt.show()