import sys
import numpy as np
import subprocess
import base64 #import datetime #import matplotlib.pyplot as plt
import os
import ctypes

mytest = 0
#=================================Input Begin===========================
if mytest: print('%Tst len(sys.argv),', len(sys.argv))
tex1 = ' *** Perimeters, areas, random inscribed polygons ***'
tex2 = '(May-202, MC)'; tex3 = '22_18:25'
print(tex1+(81-len(tex1)-len(tex2))*' '+tex2+'\n'+73*' '+tex3)

nverts, language, ntrials, iseed, klasses, ishow = sys.argv[1:]
nverts = int(nverts)
language = int(language)
ntrials = int(float(ntrials))
iseed = int(iseed)
klasses = int(klasses)
ishow = int(ishow)

if nverts < 2:
    nverts = 2; print(' %Wrn nverts has been ajusted.')
print(' {:<24}{:>12}    | (min, 2)'.format('N. of vertices,', nverts))
print(' {:<24}{:>12}    | 0, 1, 2: Python, F90, C'.format('Simulation language,', language))
print(' {:<15}{:9.3g}{:>12}    | (Random Number Generator seed)'.format('Trials, seed,', float(ntrials), iseed))

if ntrials > 1.e+7:
    ntrials = 1.e+7; print ('%Wrn Reduced to 1.e+7')
print(' {:<12}{:>24}    |'.format('Classes,', klasses))
print(' {:<24}{:>12}    |'.format('Show,', ishow))
print(' '+40*'-'+'+'+39*'-')
#=================================Input End=============================
radius = 1
peritheo = 2 * nverts * radius * np.sin(np.pi/nverts)
theta = 2 * np.pi / nverts
print (' {:<24}{:>12.4f}    | = 2\u03c0 / N; degrees,{:>6.1f}'.\
  format('\u03B8 (rad.),', theta, theta/np.pi*180))
areatheo = np.pi * radius**2 * np.sin(theta) / theta
print (' {:<24}{:>12.2f}    | 2 N R sin(\u03c0/N)'.\
  format('Max. theo. perimeter,', peritheo))
print (' {:<24}{:>12.2f}    | \u03c0 R\u00b2 sin(\u03b8)/\u03b8'.\
  format('Max. theo. area,', areatheo))

xk_vec   = np.array([0.]*klasses)
for i in range(klasses): xk_vec[i] = i * peritheo / klasses
print(' {:<16}{:>8.1f}{:>12.1f}    |'.\
  format('Abscissas in', xk_vec[0], xk_vec[-1]))
#=============================Simulation================================
#tex = ' Simulating ... '
#print(tex+(41-len(tex))*' '+'|')
x_vec = np.array([0.]*klasses)

xp_vec = np.array([0.]*klasses)
freqp_vec = np.array([0.]*klasses)
xa_vec = np.array([0.]*klasses)
freqa_vec = np.array([0.]*klasses)

if language == 0: # Python
    import ModPython
    xp_vec, freqp_vec, xa_vec, freqa_vec = ModPython.simulpolygons(nverts, iseed, ntrials, klasses)

if language == 1: # fortran
    import modRandomPolygons
    xp_vec, freqp_vec, xa_vec, freqa_vec = \
      modRandomPolygons.simulpolygons(nverts, iseed, ntrials, klasses)

if language == 2: # C
  my_c_code = ctypes.CDLL("./simulate_polygons.so")
    # Define the C function prototype
  my_c_code.simulpolygons.argtypes = [
      ctypes.c_int, ctypes.c_int, ctypes.c_int, ctypes.c_int,
      ctypes.POINTER(ctypes.c_double), ctypes.POINTER(ctypes.c_double), \
        ctypes.POINTER(ctypes.c_double), ctypes.POINTER(ctypes.c_double)]
  # Call the C function
  my_c_code.simulpolygons(
    nverts, iseed, ntrials, klasses,
    xp_vec.ctypes.data_as(ctypes.POINTER(ctypes.c_double)),
    freqp_vec.ctypes.data_as(ctypes.POINTER(ctypes.c_double)),
    xa_vec.ctypes.data_as(ctypes.POINTER(ctypes.c_double)),
    freqa_vec.ctypes.data_as(ctypes.POINTER(ctypes.c_double)))

if language > 2 or language < 0:
    print('%Err Bad language. Stop'); exit()

xp_vec = np.fromiter(xp_vec, dtype=float, count=klasses)
xa_vec = np.fromiter(xa_vec, dtype=float, count=klasses)
freqp_vec = np.fromiter(freqp_vec, dtype=float, count=klasses)
freqa_vec = np.fromiter(freqa_vec, dtype=float, count=klasses)

ypmax = np.max(freqp_vec);            yamax = np.max(freqa_vec)
xpmax = xp_vec[np.argmax(freqp_vec)]; xamax = xa_vec[np.argmax(freqa_vec)]
h = ' {:<12}{:>12.2f}{:>12.2f}    |'
print(h.format('max P at', xpmax, ypmax))
print(h.format('max A at', xamax, yamax))

#print('AQUI 3')
if nverts==2: freqa_vec = np.array([0.]*klasses)
#print (' {:<12}{:>12.4f}{:>12.4f}    | theo., {:<12.4f}{:<12.4f}'.\
#  format('Max. P, A,', perimax, areamax, peritheo, areatheo))

def averagestdev(x_vec, y_vec):
    klas = np.size(x_vec); delta = (x_vec[-1] - x_vec[0]) / klas
    ave = ((x_vec[0]*y_vec[0] + x_vec[-1]*y_vec[-1]) / 2 + \
      np.sum(x_vec[1:-2]*y_vec[1:-2])) * delta
#   ave = ((y_vec[0] + y_vec[-1]) / 2 + np.sum(y_vec[1:-2])) * delta
#   ave = ave / (x_vec[-1] - x_vec[0])
    std = (((x_vec[0]-ave)**2 * y_vec[0] + \
      (x_vec[-1]-ave)**2 * y_vec[-1]) / 2 + \
      np.sum((x_vec[1:-2]-ave)**2 * y_vec[1:-2])) * delta
    return ave, std

pave, pstd = averagestdev(xp_vec, freqp_vec)
print(' {:<12}{:>12.4g}{:>12.4g}    | Perimeter'.\
  format('P: ave std,', pave, pstd))
aave, astd = averagestdev(xa_vec, freqa_vec)
print(' {:<12}{:>12.4g}{:>12.4g}    | Area'.\
  format('A: ave std,', aave, astd))
buffer = "P ave std: %5.3g %5.3g A ave std: %5.3g %5.3g" % \
  (pave, pstd, aave, astd)

tex = ' To plot:'; print(tex+(41-len(tex))*' '+'|')
# import scipy
# print(' (a) scipy version:', scipy.__version__, 'end') # with import scipy
# from importlib.metadata import version # not on web
# print(' (b) scipy version:', version('scipy'), 'end')

basedir = '/tmp/tx'+str(os.getpid())
tmpfilexy = basedir+'tmpfxy'
filexy = open(tmpfilexy, 'w') # define 'filexy'

for i in range (klasses):
    print(xp_vec[i], freqp_vec[i], xa_vec[i], freqa_vec[i], file=filexy)

if mytest: print('%Tst To infile ...')
infile = basedir+'tmp2.dat'
imgfile = basedir+'tmp2.png'
open(imgfile, 'w')

filexy.seek(0)
my_loadstr = \
f'''wid = 620; set term png font 'Times, 14' size wid, wid/1.618
set output '{imgfile}'
set title 'Perim., Area, nv = {nverts}' offset 0, -0.35
set key center right Right
set timestamp font 'Helvetica-Bold, 9'
set xlabel '{{/Times-Italic p}}' offset +8., 0.75
set x2label '{{/Times-Italic a}}' offset +8., -0.4
set xtics nomirror; set x2tics nomirror
set ylabel 'pdf_P' rotate by 90 offset 1.5, 0
set y2label 'pdf_A' rotate by 90 offset -2.5, 0
set ytics nomirror; set y2tics nomirror
set label '{buffer}' at character 12, character 7
set style line 12 lc rgb 'red' lw 3 pt 21 ps 2
set style line 13 lc rgb 'black' lw 1 pt 19 ps 1.2
plot '{tmpfilexy}' usi 1:2 axis x1y1 title 'pdf_P' w lines ls 12, \
 '' usi 3:4                axis x2y2 title 'pdf_A' w lines ls 13
'''
if mytest:
   print('%Tst To print myloadstr ...'); print(my_loadstr)

with open(infile, 'w') as pli:
    lines = my_loadstr.splitlines(True)
    pli.writelines(lines)
    pli.flush()
    subprocess.call(['/usr/bin/gnuplot', infile], shell=False)
###s = os.system('/usr/bin/wc '+imgfile)

with open(imgfile, 'rb') as image_file:
    encoded = base64.b64encode(image_file.read())
encoded = str(encoded)[2:-1]
os.remove(infile); os.remove(imgfile); os.remove(tmpfilexy)

print('<center><!--PLOT:-->')
print('<img alt="Plot" src="data:image/png; base64, '+encoded+'" />')
print('</center>')

if ishow:
    print(' Show_i___________x__________yy...')
    for i in range(klasses):
        print(' {:6}{:12.6g}{:12.6g}{:12.6g}{:12.6g}'.\
          format(i, xp_vec[i], freqp_vec[i], xa_vec[i], freqa_vec[i]))

print(' End of program')