""" From "COMPUTATIONAL PHYSICS" & "COMPUTER PROBLEMS in PHYSICS"
    by RH Landau, MJ Paez, and CC Bordeianu (deceased)
    Copyright R Landau, Oregon State Unv, MJ Paez, Univ Antioquia, 
    C Bordeianu, Univ Bucharest, 2018. 
    Please respect copyright & acknowledge our work."""

# CatFrictionVis.py: Visual, waves on catenary with friction

import numpy as np

dt = 0.0001
dx = 0.01
T = 1
rho = 0.1
maxtime = 100
kappa = 30
D = T / (rho * 9.8)

x = np.zeros((512, 3), float)
q = open('CatFriction.dat', 'w')
rr = open('CatFunct.dat', 'w+t')

for i in range(0, 101):
    x[i][0] = -0.08 * np.sin(np.pi * i * dx)  # IC

for i in range(1, 100):              # First step
    x[i][1] = (dt * (T / rho) * ((x[i + 1][0] - x[i][0])
                                 / dx * (np.exp((i - 50) * dx / D)
                                         - np.exp(-(i - 50) * dx / D)) / D + (np.exp((i - 50) * dx / D)
                                                                              + np.exp(-(i - 50) * dx / D)) * (x[i + 1][0] + x[i - 1][0]
                                                                                                               - 2.0 * x[i][0]) / (pow(dx, 2)))
               - 2 * kappa * x[i][0] + 2 * x[i][0] / dt) / (2 * kappa + (2 / dt))

for k in range(0, 300):             # Other steps
    for i in range(1, 100):
        x[i][2] = (dt * (T / rho) * ((x[i + 1][1] - x[i][1])
                                     / dx * (np.exp((i - 50) * dx / D) - np.exp(-(i - 50) * dx / D)) / D
                                     + (np.exp((i - 50) * dx / D) + np.exp(-(i - 50) * dx / D)) *
                                     (x[i + 1][1] + x[i - 1][1] - 2.0 * x[i][1]) / (pow(dx, 2)))
                   - 2 * kappa * x[i][1] - (-2 * x[i][1] + x[i][0]) / dt) / (2 * kappa + (1 / dt))
    for i in range(1, 101):
        x[i][0] = x[i][1]
        x[i][1] = x[i][2]
    if (k % 4 == 0 or k == 0):
        for i in range(0, 100):
            a1 = np.exp((i - 50.) * dx / D)
            a2 = np.exp(-(i - 50.) * dx / D)
            rr.write("%7.3f" % (D * (a1 + a2)))
            rr.write("\n")
            q.write('%7.3f' % (x[i, 2]))
            q.write("\n")
        q.write("\n")
        rr.write("\n")
rr.closed
q.closed
print("Data stored in CatFrict.dat and CatFunct.dat")
