######################################
# MAIN PROGRAM
# creates compilable FORTRAN code
# based on parameter.in and sis_src.f
# compiles and executes it in parallel
######################################

import fileinput
import numpy as np
import pandas as pd
import os
import subprocess
from multiprocessing.dummy import Pool as ThreadPool

######################################
# executes bash commands and 
# communicates output (None,None)
######################################
def bash_command(cmd):
    p=subprocess.Popen(['/bin/bash', '-c', cmd],stderr=subprocess.PIPE)
    print p.communicate(),cmd

######################################
# converts numeric values to fortran's
# double precision format
######################################
def givedouble(var):
    return ("%.1e"%var).replace("e","d")

#####################################
# replaces text inline for a specific
# search string (sext->sfor) in ffile
#####################################
def myreplace(sline,sext,sfor,ffile,found):
    if(sline.find(sext)>=0):
        ffile.write(sline.replace(sext, sfor)),
        return found+1
    else:
        return found+0

######################################
# creates compilable fortran code from
# the source code in sis_src.f and our
# input parameters passed by the main
# program but generated by input.py
######################################
def createdotf(pars):
    identifier,fname,fstring,md,ac,ms,metal,tgas,gama,cmigI,Apert,Fpert,Tfin,Verb = pars

# calculates toomre q (stability criterion is q>1)
    amin=ac*((7/4.-gama)/(2.-gama))**(1./(2-gama))
    sigmag=(2.-gama)*md*2.e33/(2.*np.pi*(ac*1.5e13)**2.)
    qc=1.24e5*(amin**(gama-7./4.))*(ac**(-gama))*(ms)*np.exp((amin/ac)**(2.-gama))/sigmag
# convert input to fortran formatting
    cmigI=givedouble(cmigI)
    Fpert=givedouble(Fpert)
    md=givedouble(md)
    Tfin=givedouble(Tfin)
    Verb=".%s."%str(Verb).upper()
# checks for disk stability, it shouldn't be a problem if parameter.in is obtained from input.py   
    if(qc >= 1.):
        f = open(fstring+fname+'%s.f'%identifier, 'w')
        for line in fileinput.input("sis_src.f", inplace=False):
            found=0
    # output filenames (do not change, rescortos is called by pastefiles.sh):
    # resultados - log time evolution (only when verbose=true)
    # rescortos - default results
    # gigantes - info for giant planet formation
            found=myreplace(line,"$ArchivoSalida","\""+fstring+"resultados%s.txt\""%identifier,f,found)
            found=myreplace(line,"$SalidaCorta","\""+fstring+"rescortos%s.txt\""%identifier,f,found)
            found=myreplace(line,"$Gigantes","\""+fstring+"gigantes%s.txt\""%identifier,f,found)
    # simulation parameters:
    # disk mass
    # characteristic radius (disk)
    # stellar mass
    # toomre q parameter
    # metallicity
    # gas dissipation timescale
            found=myreplace(line,"$Emed",str(md),f,found)
            found=myreplace(line,"$Rcar",str(ac),f,found)
            found=myreplace(line,"$Emestar",str(ms),f,found)
            found=myreplace(line,"$Qest",str(qc),f,found)    
            found=myreplace(line,"$Emetal",str(metal),f,found) 
            found=myreplace(line,"$Taugas",str(tgas),f,found) 
    # sort-of-fixed parameters:
    # power law profile for surface density
    # type I migration parameter
            found=myreplace(line,"$Gama",str(gama),f,found)
            found=myreplace(line,"$ConstMigI",str(cmigI),f,found)
    # perturbation parameters:
    # amplitude
    # length scale
            found=myreplace(line,"$Apert",str(Apert),f,found)
            found=myreplace(line,"$Fpert",str(Fpert),f,found)
    # simulation timescale
            found=myreplace(line,"$Tfin",str(Tfin),f,found)
    # true for verbose output (log-time evolution)        
            found=myreplace(line,"$Verb",str(Verb),f,found)
    # add identifier to output
            found=myreplace(line,"$Identifier",str(identifier),f,found)
            if(found==0):
                f.write(line)
        f.close()
# writes simulation parameters into info file    
    finfo = open(fstring+'info%s.txt'%identifier, 'w')
    firstline="%s, %s, %s, %s, %s, %s, %s, %s, %s, %s, %s, %s\n\n" % (identifier,md,ac,ms,metal,tgas,gama,cmigI,Apert,Fpert,Tfin,Verb)
    finfo.write(firstline)
    finfo.write("Identifier...............................%s\n"%identifier)
    finfo.write("Disk Mass (M_sun)........................%s\n"%md)
    finfo.write("Critical Radius (AU).....................%s\n"%ac)
    finfo.write("Stellar Mass (M_sun).....................%s\n"%ms)
    finfo.write("Metallicity (log Fe/H)_*/(log Fe/H)_sun..%s\n"%metal)
    finfo.write("Gas depletion timescale (yr).............%s\n"%tgas)
    finfo.write("Gamma (density profile power law)........%s\n"%gama)
    finfo.write("Type I Migration Coefficient.............%s\n"%cmigI)
    finfo.write("Perturbation Amplitude...................%s\n"%Apert)
    finfo.write("Perturbation length scale................%s\n"%Fpert)
    finfo.write("Simulation timescale (yr)................%s\n"%Tfin)
    finfo.write("Verbose results?.........................%s\n"%Verb)
    finfo.write("Toomre Q.................................%s\n"%qc)
    finfo.close()
    
#   reads from parameter.in, which should be generated from input.py
df=pd.read_csv('parameters.in')

cstring1 = 'gfortran -o' # compiler+flags
fname = 'qsistemas' # file names
cstring3 = '.f; ./' # extensions (do not change)
fstring  = 'results/' # folder for .f sources, executables and results (do not change, req. by pastefiles.sh)

# creates work folder
if not os.path.exists(fstring):
    os.makedirs(fstring)

nsys=len(df) # number of simulations
myfill=int(np.log10(nsys)+1) # with zfill, generate identifier for each simulation
lidr=[]
inst=[]
for icoun in range(nsys):
    idr=str(icoun).zfill(myfill)
    pars=tuple([idr]+[fname]+[fstring]+df.iloc[icoun].values.tolist()) # parameters for createdotf function
    createdotf(pars)
    pname=fname+idr    
    inst.append(cstring1+2*(' '+fstring+pname)+cstring3+fstring+pname) # generate instructions for bash_command function
    
# create header for csv (do not change filename, req. by pastefiles.sh)
hfile=open(fstring+'header.txt', 'w')
hfile.write('ident,it,t,a(i),emegas(i),emepla(i)/emet,rplanet(i)/radtie,emestar,rc,qest,sigmag_0,emed,gama,apert,fpert,constmigI,emetal,taugas\n') 
hfile.close()

# parallel parsing of instructions (compile+execute from synthetic .f sources) to pool of workers
pool = ThreadPool(5)
results=pool.map(bash_command,inst)
pool.close()
pool.join()

# open parameters in sis_src.f:
# lsnr : number of log time snapshots for verbose output (500 yields about 2-5 MB)
# dt: timestep in years (default to 100)

# todo:
# perturbation only works for the dust - maybe add gas perturbation?
# gas envelope is automatically dissipated after a collision (look for eme1 and eme2) - modulate better?
# make a verbose option for giant planets
# resultados files go up tfin in parameters.in, but if verbose=false this is overkill (curb tfin to 101% in src?)
# auto check empty files and communicate in paste bash script

# bugs:
# not so much a bug, but old input is not trustworthy, as the q condition was not calculated properly
# i for kk in planet dissapearance when w > wcrit, caused whole systems to disappear due to bad indexing
# massive removal of runtime errors due to segfaults (april '16)

# notes:
## an empty rescortos file (quite possibly) happens because of migration + collisions (loss of env) + dissasembly due to high L
#  when this happens, check gigantes file for giant planet formation + migration
## sis_src outputs at 2e7 yr for rescortos by default

