# derived from database_write.py 
# writes all the 4 database files suite in one run.
# first used to add radiation peaking data

# V02 with treatment of m2 cases
# writes a datafile for given:
# nmode, NBIthreshold, radiation threshold, MODE threshold
# one line for each pulse in list0 with NBI threshold exceeded.
# The radiation threshold condition is: 
# Prad > 0.7*Ptot & t > NBIstart + 0.3s (to avoid initial NBI transients)
# it is stupidly hidden on the powers function of databasefunctions.py
# Mode recognized if at least 12 points (100 ms) are above threshold
#npulse   # pulse
#B0       # B0 T
#Ip       # Ip MA  mean on Pt1-Pt2 interval
#Pt1 Pt2  # NBI on an off times (threshold)
#t_rad    # time of radiation crossing a thershold fraction of total power
#P        # max NBI power MW                                                           <<< Pmax
#...      # relative power variation in shorter interval after mode on (def. 0.2 s)   <<<<<<<<  ADD <<< dP1
#...      #      "                   in longer interval after onset (def. 0.4 s)      <<<<<<<<  ADD <<< dP2

#q95      # q95, 
#tril     # lower triangularity
#triu     # upper triang
#bNm      # max betaN                                                                 <<< bNmax
#bNt      # time of max betaN                                                         <<< tbNmax
#bNon     # betaN at mode onset
#...      # beta var. at the end of shorter interval after mode on (def. 0.2 s)      <<<<<<<<  ADD <<< dbN1
#...      # beta var. at the end of longer interval after onset (def. 0.4 s)         <<<<<<<<  ADD <<< dbN2

#Mmax     # Mode maximum amplitude before t_rad  (0 if no mode)
           #Also available max amplitude after t_rad (mode in tail)
#Mdur     # Mode duration to last time above threshold   (0 if no mode)          
#ONt      # Mode onset time   (0 if no mode)                                          <<< Mtime
#ONtype   # rapidity of intial growh  (10 if no mode)         NOT WELL DEFINED        <<< Mtype
#ONa      # onset amplitude: mean on 12 pts after threshold crossing  (0 if no mode)  <<< Mamp0
           #the condition for mode existence is at least 12 pts above threshold. POSSIBLE FALSE ALARM
#...      # average amplitude in shorter interval after onset (def. 0.2 s)      <<<<<<<<  ADD <<< Mamp1
#...      # average amplitude in longer interval after onset (def. 0.4 s)       <<<<<<<<  ADD <<< Mamp2
 
#kg1l     # line average density from interferomery
#hrtx     # line average density from HRTS

import numpy as np
import scipy as sp
import matplotlib.pyplot as plt
import ppfeg
import databasefunctions as df



# IMPROVE USE AND DESCRIPTION OF   dt1 and dt2
# WRITE INPUTS IN A HEADER
# TRY FIXED FORMATS

def producefile(pulselist, outfile, NBIthresh, dtthresh, dt1, dt2, dtx, maxfrad, nmode, m2, MODEthresh):
	# file header with options
	headformats = ', NBIthresh={0}MW, dtthresh={1}s, dt1={2}s, dt2={3}s, dtx={4}s, maxfrad={5}, Modethresh={6}G'
	firstline = pulselist + headformats.format(1.e-6*NBIthresh, dtthresh,\
		dt1, dt2, dtx, maxfrad, 1.e4*MODEthresh)
	# n1m2 exceptions for hybrids
	m2sole = [73031, 73034, 73308, 73310, 73725, 73729, 73752, 73754, 73755, 73756, 74190, 75389]
	m2sole.extend([83413, 83417, 83521, 84675])
	m2after = {73749:48.2, 73751:47.58, 74184:47.81, 74185:47.13, 74186:50.13, 74195:50.3, 74196:47.35, 75584:46.13}
	m2after.update({83528:46.6, 83564:45.99, 85204:45.41, 87014:45.25})
	
	names = 'npulse B0    Ip   '                          # pulse numner, B0 (T), Ip (MA)
	formats = '{0} {1:5.2f} {2:5.2f} '
	names = names + ' Pt1     Pt2     t_rad  Pmax '            # NBI on and off t, radiation time, max NBI (MW)
	formats = formats + '{3:7.3f} {4:7.3f} {5:7.3f} {6:6.3f} '
	names = names + '   Pon     dP1    dP2 '          # NBI power at onset and rel. variation in 0.2, 0.4 s
	formats = formats + '{7:7.3f} {8:6.3f} {9:6.3f} '
	names = names + '   kg1l  hrtx '                     # line average density from interf. and HRTS
	formats = formats + '{10:5.2f} {11:5.2f} '
	names = names + ' q95   tril  triu  bNmax tbNmax '    # q95, triangs, max betaN, tmax
	formats = formats + '{12:5.2f} {13:5.2f} {14:5.2f} {15:5.2f} {16:7.3f} '
	names = names + ' bNon   dbN1  dbN2 '                # betaN @ onset and variation 0.2, 0.4 s after
	formats = formats + '{17:6.3f} {18:5.2f} {19:5.2f} '
	names = names + ' Hon    dH1   dH2 '     # H98 at onset and variation 0.2, 0.4 s after
	formats = formats + '{20:6.3f} {21:5.2f} {22:5.2f} '
	names = names + '  Mmax   Mdur  Mtime  Mtype '         # Mode: max, duration, on_time, on_type
	formats = formats + '{23:6.3f} {24:6.3f} {25:6.3f} {26:4.1f} '
	names = names + 'Mamp0  Mamp1  Mamp2 '             # Mean amplitude 0.1, 0.2, 0.3 s after onset
	formats = formats + '{27:6.3f} {28:6.3f} {29:6.3f} '
	names = names + '  BOLOon   dBOLO1  dBOLO2 '             # Core radiation
	formats = formats + '{30:7.0f} {31:8.0f} {32:8.0f} '
	names = names + 'xpeakon dxpeak1 dxpeak2 '             # Soft-x peaking
	formats = formats + '{33:5.2f} {34:8.2f} {35:7.2f} '
	names = names + 'Pt1s   Pt2s   Pmaxs '             # power data with smoothing
	formats = formats + '{36:7.3f} {37:7.3f} {38:6.3f} '
	names = names + '   kg1v  li3m '                     # line int. density and internal inductance
	formats = formats + '{39:5.2f} {40:5.2f} '
	formats = formats + '\n'
	
	fo = open(outfile, 'w')
	fo.write(firstline + '\n')
	fo.write(names + '\n')
	
	with open(pulselist) as f:
		for line in f:
			npulse = int(line)
			print npulse
			dda = 'QMHD'
			if npulse < 80000: dda = 'CMHD'
			
			#Pt1, Pt2, t_rad, P, timebolo, nbi, ptot, bolotopi = df.powers(npulse, NBIthresh, maxfrad=maxfrad)
			Pt1, Pt2, t_rad, P, timebolo, nbi, ptot, bolotopi, Pt1s, Pt2s, Ps = \
				df.powers5(npulse, NBIthresh, dtthresh=dtthresh, maxfrad=maxfrad)
			if Pt1 < 41.0:
				continue
	
			Ip = df.IMA(npulse, Pt1, t_rad)
			
			B0 = df.btor(npulse, Pt1, t_rad)
			
			if ( nmode == 1 & ~m2 ):           # n=1 m=1 mode
				if ( npulse in m2sole ):
					Mtail, Mmax, Mdur, Mtime, Mtype, M0, M1, M2 = \
						df.newMODE(npulse, 10, MODEthresh, Pt1, t_rad, dda=dda, dt1=dt1, dt2=dt2)
				elif ( npulse in m2after ):
					Mtail, Mmax, Mdur, Mtime, Mtype, M0, M1, M2 = \
						df.newMODE(npulse, nmode, MODEthresh, Pt1, m2after[npulse], dda=dda, dt1=dt1, dt2=dt2)
						#df.newMODE(npulse, nmode, MODEthresh, Pt1, t_rad, dda=dda, dt1=dt1, dt2=dt2)
					#Mdur = m2after[npulse] - Mtime
				else:
					Mtail, Mmax, Mdur, Mtime, Mtype, M0, M1, M2 = \
						df.newMODE(npulse, nmode, MODEthresh, Pt1, t_rad, dda=dda, dt1=dt1, dt2=dt2)
						
			if ( nmode == 1 & m2 ):            # n=1 m=2 mode
				if ( npulse in m2sole ):
					Mtail, Mmax, Mdur, Mtime, Mtype, M0, M1, M2 = \
						df.newMODE(npulse, nmode, MODEthresh, Pt1, t_rad, dda=dda, dt1=dt1, dt2=dt2)
				elif ( npulse in m2after ):
					t_start = m2after[npulse]
					Mtail, Mmax, Mdur, Mtime, Mtype, M0, M1, M2 = \
						df.newMODE(npulse, nmode, MODEthresh, t_start, t_rad, dda=dda, dt1=dt1, dt2=dt2)
				else:
					Mtail, Mmax, Mdur, Mtime, Mtype, M0, M1, M2 = \
						df.newMODE(npulse, 10, MODEthresh, Pt1, t_rad, dda=dda, dt1=dt1, dt2=dt2)
				
			if nmode > 1:                      # other modes
				Mtail, Mmax, Mdur, Mtime, Mtype, M0, M1, M2 = \
					df.newMODE(npulse, nmode, MODEthresh, Pt1, t_rad, dda=dda, dt1=dt1, dt2=dt2)
			
			bNmax, tbNmax, bNon, dbN1, dbN2 = df.beta(npulse, Mtime, dt1=dt1, dt2=dt2) 
			Hon, dH1, dH2 = df.H98(npulse, Mtime, tbNmax, dt1=dt1, dt2=dt2)
			BOLOon, dBOLO1, dBOLO2, xpeakon, dxpeak1, dxpeak2 = \
				df.radiation(npulse, Mtime, tbNmax, dt1=dt1, dt2=dt2, dtx=dtx)
			Pon, dP1, dP2 = df.pvar(npulse, Mtime, timebolo, nbi, dt1=dt1, dt2=dt2)
			
			q95, tril, triu = df.q_95(npulse, Mtime, tbNmax)
			
			kg1l, hrtx = df.dens(npulse, Mtime, tbNmax)
			kg1v, li3m = df.dens5(npulse, Mtime, tbNmax)
			#print Pt1s, Pt2s, Ps, kg1v, li3m
			
			fo.write(formats.format(npulse, B0, Ip, float(Pt1), float(Pt2), float(t_rad), 
				P, Pon, dP1, dP2, 1.e-19*kg1l, 1.e-19*hrtx, float(q95), float(tril), float(triu), 
				bNmax, tbNmax, bNon, dbN1, dbN2, Hon, dH1, dH2, 
				1.e4*Mmax, Mdur, Mtime, Mtype, 1.e4*M0, 1.e4*M1, 1.e4*M2, 
				BOLOon, dBOLO1, dBOLO2, xpeakon, dxpeak1, dxpeak2, 
				float(Pt1s), float(Pt2s), Ps, 1.e-19*kg1v, float(li3m)))

	fo.close()	

# General options
dt1 = -0.3             # negative to detrend
dt2 = 0.3
dtx = 0.05             # interval for soft-x peaking variation
maxfrad = 0.7          # maximum accepted radiation fraction
#NBIthresh = 4.5e6
MODEthresh = [0.0, 1.e-4, 3.e-5, 1.e-5]

#New thresholds
NBIthresh = 8.0e6
dtthresh = 0.45


# files with all data at max betaN (by large MODEthresh=100.0)
nmode = 1
letter = 'y'
#pulselist = 'list_test' 
pulselist = 'list_hybrids_ILW_sel'
m2 = False; outfile = 'n1_' + letter + '_ILW'         # n=1, m=1
#pulselist = 'list_hybrids_CW_sel2'
#m2 = False; outfile = 'n1_' + letter + '_CW'         # n=1, m=1
producefile(pulselist, outfile, NBIthresh, dtthresh, dt1, dt2, dtx, maxfrad, nmode, m2, 100.0)

"""

letter = 'f'
# ILW hybrids
#pulselist = 'list_test' 
pulselist = 'list_hybrids_ILW_sel'
nmode = 1
m2 = False; outfile = 'n1_' + letter + '_ILW'         # n=1, m=1
producefile(pulselist, outfile, NBIthresh, dtthresh, dt1, dt2, dtx, maxfrad, nmode, m2, MODEthresh[nmode])
m2 = True; outfile = 'n1m2_' + letter + '_ILW'        # n=1, m=2
producefile(pulselist, outfile, NBIthresh, dtthresh, dt1, dt2, dtx, maxfrad, nmode, m2, MODEthresh[nmode])
nmode = 2
outfile = 'n2_' + letter + '_ILW'                     # n=2
producefile(pulselist, outfile, NBIthresh, dtthresh, dt1, dt2, dtx, maxfrad, nmode, m2, MODEthresh[nmode])
nmode = 3
outfile = 'n3_' + letter + '_ILW'                     # n=3
producefile(pulselist, outfile, NBIthresh, dtthresh, dt1, dt2, dtx, maxfrad, nmode, m2, MODEthresh[nmode])

# CW hybrids
# NO DATA FROM SOFT-X
pulselist = 'list_hybrids_CW_sel2'
nmode = 1
m2 = False; outfile = 'n1_' + letter + '_CW'         # n=1, m=1
producefile(pulselist, outfile, NBIthresh, dtthresh, dt1, dt2, dtx, maxfrad, nmode, m2, MODEthresh[nmode])
m2 = True; outfile = 'n1m2_' + letter + '_CW'        # n=1, m=2
producefile(pulselist, outfile, NBIthresh, dtthresh, dt1, dt2, dtx, maxfrad, nmode, m2, MODEthresh[nmode])
nmode = 2
outfile = 'n2_' + letter + '_CW'                     # n=2
producefile(pulselist, outfile, NBIthresh, dtthresh, dt1, dt2, dtx, maxfrad, nmode, m2, MODEthresh[nmode])
nmode = 3
outfile = 'n3_' + letter + '_CW'                     # n=3
producefile(pulselist, outfile, NBIthresh, dtthresh, dt1, dt2, dtx, maxfrad, nmode, m2, MODEthresh[nmode])
"""

"""
# ILW baseline
pulselist = list_baseline_ILW_sel
nmode = 1
m2 = False; outfile = 'n1_' + letter + '_baseline_ILW'         # n=1, m=1
producefile(pulselist, outfile, NBIthresh, dt1, dt2, dtx, maxfrad, nmode, m2, MODEthresh[nmode])
m2 = True; outfile = 'n1m2_' + letter + '_baseline_ILW'        # n=1, m=2
producefile(pulselist, outfile, NBIthresh, dt1, dt2, dtx, maxfrad, nmode, m2, MODEthresh[nmode])
nmode = 2
outfile = 'n2_' + letter + '_baseline_ILW'                     # n=2
producefile(pulselist, outfile, NBIthresh, dt1, dt2, dtx, maxfrad, nmode, m2, MODEthresh[nmode])
nmode = 3
outfile = 'n3_' + letter + '_baseline_ILW'                     # n=3
producefile(pulselist, outfile, NBIthresh, dt1, dt2, dtx, maxfrad, nmode, m2, MODEthresh[nmode])
"""
