#! /usr/bin/env python

# needs python, gcc, gnuplot, latex, dvips, ps2pdf, ps2eps

import math
import os

def gnuplot_set_style_line(line_id = '1', linetype = '1', linecolor = '7', linewidth = '1', pointtype = '7', pointsize = '1'):
    """Set a gnuplot line style."""
    return ('set style line ' + line_id   + ' lt ' + linetype
                     + ' lc ' + linecolor + ' lw ' + linewidth
                     + ' pt ' + pointtype + ' ps ' + pointsize + '\n')

def turn_epslatex_to_eps_plain(use_euler, filename, figname):
    tmp_tex = open('turn_epslatex_to_eps_tmp.tex','w')
    tmp_tex.write('\\documentclass[11pt]{article}\n'
                + '\\usepackage{amsmath, amsfonts, amssymb}\n'
                + '\\usepackage{color}\n'
                + '\\usepackage{graphicx}\n')
    if use_euler == 1:
        tmp_tex.write('\\usepackage[mathcal, mathbf]{euler}\n')
    tmp_tex.write('\\pagestyle{empty}\n'
                + '\\begin{document}\n'
                + '\\begin{figure}\n'
                + '\\begin{center}\n'
                + '\\input{' + filename + '}\n'
                + '\\end{center}\n'
                + '\\end{figure}\n'
                + '\\end{document}\n')
    tmp_tex.close()
    os.system('latex turn_epslatex_to_eps_tmp.tex')
    os.system('dvips turn_epslatex_to_eps_tmp.dvi')
    os.system('ps2eps -f turn_epslatex_to_eps_tmp.ps')
    os.system('mv turn_epslatex_to_eps_tmp.eps ' + figname + '.eps')
    os.system('rm turn_epslatex_to_eps_tmp.ps '
               + 'turn_epslatex_to_eps_tmp.dvi '
               + 'turn_epslatex_to_eps_tmp.tex '
               + 'turn_epslatex_to_eps_tmp.aux '
               + 'turn_epslatex_to_eps_tmp.log')
    return None

def plot_epslatex_plain(plot_command, linewidth, sizex, sizey, use_euler, figname):
    out = open(figname + '_gnuplot_instructions.txt','w')
    out.write('set term epslatex color lw {0} size {1},{2} font ",9.1"\n'.format(linewidth, sizex, sizey))
    out.write('set key spacing 1.5\n')
    out.write('set output "' + figname + '_tmp.tex"\n')
    out.write(plot_command)
    out.write('set output\n')
    out.close()
    os.system('rm ' + figname + '.eps')
    gnuplot_result = os.system('gnuplot ' + figname + '_gnuplot_instructions.txt')
    if gnuplot_result == 0:
        os.system('rm ' + figname + '_gnuplot_instructions.txt')
        turn_epslatex_to_eps_plain(use_euler, figname + '_tmp', figname)
        os.system('rm ' + figname + '_tmp.eps ' + figname + '_tmp.tex')
    else:
        print '********\ngnuplot error for figure ' + figname + '\n********\n'
    return gnuplot_result

class bifurcation_computer:
    def __init__(self, name, formula, rmin, rmax, xmin, xmax):
        self.name = name
        self.formula = formula
        self.rmax = rmax
        self.rmin = rmin
        self.xmax = xmax
        self.xmin = xmin
        return None
    def get_cycles(self, random_seed, parameter_resolution,
                   number_of_points, number_of_iterations, iterations_to_keep):
        print self.formula
        src = open(self.name + '.c', 'w')
        src.write('#include <stdio.h>\n'
                + '#include <stdlib.h>\n'
                + '#include <math.h>\n\n'
                + 'int main()\n{\n'
                + '    int i, j, k;\n'
                + '    double x, r;\n'
                + '    srand48({0});\n'.format(random_seed)
                + '    for (k=1; k<={0}; k++)\n'.format(parameter_resolution)
                + '    {\n'
                + '        r = {0} + {1}*(double)(k)/{2};\n'.format(self.rmin, self.rmax - self.rmin, parameter_resolution)
                + '        for (i=0; i<{0}; i++)\n'.format(number_of_points)
                + '        {\n'
                + '            x = {0} + drand48()*{1};\n'.format(self.xmin, self.xmax - self.xmin)
                + '            for (j=0; j<({0}-{1}); j++)\n'.format(number_of_iterations, iterations_to_keep)
                + '                x = ' + self.formula + ';\n'
                + '            for (; j<{0}; j++)\n'.format(number_of_iterations)
                + '            {\n'
                + '                x = ' + self.formula + ';\n'
                + '                printf("%g\\t%g\\n", r, x);\n'
                + '            }\n'
                + '        }\n'
                + '    }\n'
                + '    return 0;\n}\n')
        src.close()
        os.system('gcc -O3 -mtune=native -lm ' + self.name + '.c -o ' + self.name)
        os.system('./' + self.name + ' > ' + self.name + '.out_first')
        os.system('sort -u ' + self.name + '.out_first > ' + self.name + '.out')
        os.system('rm ' + self.name + '.c ' + self.name + '.out_first ' + self.name)
        return None

def main():
    twopi = 4*math.acos(.0)
    bc = bifurcation_computer('sinmap', 'r*sin(x)', -twopi, twopi, -twopi, twopi)
    bc.get_cycles(31, 2000, 128, 5000, 2)
    misc_command = ('set key bottom center\n'
                  + 'set samples 400\n'
                  + 'set xlabel "$r$"\n'
                  + 'set ylabel "$x$"\n'
                  + 'set key width -15\n'
                  + gnuplot_set_style_line(line_id = '1', linetype = '1', linecolor = '1')
                  + gnuplot_set_style_line(line_id = '2', linetype = '1', linecolor = '2')
                  + gnuplot_set_style_line(line_id = '3', linetype = '2', linecolor = '1')
                  + gnuplot_set_style_line(line_id = '4', linetype = '2', linecolor = '2'))
    plot_command = ('plot [{0}:{1}][{0}:{1}] "sinmap.out" notitle w d lc 7,'.format(-twopi, twopi)
                     + 'x*sin(x) title "$r \\\\sin\\\\left(r\\\\right)$" ls 1,'
                     + '-x*sin(x) title "$- r \\\\sin\\\\left(r\\\\right)$" ls 2,'
                     + 'x title "$r$" ls 3, -x title "$-r$" ls 4\n')
    plot_epslatex_plain(misc_command + plot_command, 1.0, 6, 3, 1, './sinmap')
    misc_command = ('set key bottom center\n'
                  + 'set key width -15\n'
                  + 'set xlabel "$x$"\n'
                  + 'set ylabel "$y$"\n'
                  + 'g(x,y) = x*sin(y)\n'
                  + 'set samples 1000\n')
    plot_command = ('plot [.0:{0}] "sinmap.out" notitle w d lc 7,'.format(twopi*.75)
                  + 'g(x, x) notitle w l ls 1, -g(x, x) notitle w l ls 1, g(x, x) notitle w l ls 1,'
                  + 'g(x, g(x, x)) notitle w l ls 1, -g(x, g(x, x)) notitle w l ls 1,'
                  + 'g(x, g(x, g(x, x))) notitle w l ls 1, -g(x, g(x, g(x, x))) notitle w l ls 1,'
                  + 'g(x, g(x, g(x, g(x, x)))) notitle w l ls 1, -g(x, g(x, g(x, g(x, x)))) notitle w l ls 1\n')
    plot_epslatex_plain(misc_command + plot_command, 1.0, 6, 3, 1, './envelopes')
    bc = bifurcation_computer('logmap', 'r*x*(1.0-x)', .0, 4., .0, 1.)
    bc.get_cycles(11, 2000, 128, 5000, 2)
    misc_command = ('set xlabel "$r$"\n'
                  + 'set ylabel "$x$"\n')
    plot_epslatex_plain(misc_command + 'plot "logmap.out" notitle w d lc 7\n', 1.0, 6, 3, 1, './logmap')
    return None

if __name__ == '__main__':
    main()