#!/usr/bin/env python3
# -*- coding: utf-8 -*-
import numpy as np
import matplotlib.pyplot as plt
import random
import sys

# Initialize
plt.close('all')

# calculs de trajectoire parabolique
def parabole_fc(xx, angled, vini, y0):
    # g = 5.8  # Acceleration due to gravity in m/s^2
    g = 9.8  # Acceleration due to gravity in m/s^2
    yy = y0 + np.tan(np.radians(angled)) * xx - (g / 2) * ((xx / (vini * np.cos(np.radians(angled))))**2)
    return yy

# Implémentation du maillage
xmax = 400
dx = 2
x = np.arange(0, xmax + dx, dx)
Nx = len(x) # xmax/dx +1
# print(Nx)
# sys.exit()
n_iterations = 1000
h = np.zeros((Nx))

# Building the topography of the cinder cone
n_plots = 50
i_plot = n_iterations // n_plots
j = 0

# Initialisation des paramètres
dh = .8
theta_min = 75 # 75
theta_max = 89  # 90
V0 = 50.0 # vitesse initiale 50 +- 10
dV0 = 10.0
hmax = 200
crit_slope = np.tan(np.radians(30)) # pente critique

# test taille particules
if dh/dx > crit_slope:
    print('particules trop grande ou dx trop faible')
    sys.exit()

# Create figures
fig = plt.figure(0, figsize=(10,8))
ax = fig.add_subplot()
ax.set_xlim(-xmax, xmax)
ax.set_ylim(0, hmax)
ax.set_title('Growth of a Cinder Cone without landsliding', fontsize=18)
ax.set_xlabel('Distance from vent (m)', fontsize=18)
ax.set_ylabel('Height of cinder cone (m)', fontsize=18)

for i in range(1, n_iterations + 1):
    v0 = random.gauss(V0, dV0)
    angle0 = random.uniform(theta_min, theta_max)

    # Calculate trajectory
    y = parabole_fc(x, angle0, v0, h[0])

    # Find where the bomb hits
    # print(np.argmax(y < h[:,i]))
    # print(np.where(y < h[:,i]))
    # sys.exit(0)
    landing_point = np.argmax(y < h)

    if 0 < landing_point < Nx - 1:
      h[landing_point] += dh

    if i % i_plot == 0:
        print('Iteration: {}'.format(i))
        colors = ['black', 'lightcoral','coral','peru','orange','yellow','lightgreen','green','cyan','darkviolet']
        ax.plot(x, y, c='red', alpha=0.1)
        print(j, j%len(colors))
        ax.plot(x, h, color = colors[j%len(colors)], alpha= 0.4)
        ax.plot(-x, h, color=colors[j%len(colors)], alpha=0.4)
        j = j+1
        plt.pause(0.05)
        plt.draw()

plt.show()
# Clean up
plt.close()