# -*- coding: utf-8 -*-
"""
Created on Sat Mar 14 18:52:48 2015

@author: Ondřej Ficker

Dose estimate script for Golem discharges. The full sample frequency (250 MHz) data of
HXR detector + NI basic DAS must be used (otherwise change the settings for
time axis), the data are loaded from binary file named by shot number
w/o extension.

This script will be connected to the GDB

"""
#import sys
import numpy as np
import matplotlib.pyplot as plt
import scipy.integrate as scint



sh=18979 #setting shot number
U=525    #setting operating volteg of HXR detector power source
CsE=662  #setting characteristic peak energy for callibration (in keVs)
mscint=0.38 #setting mass of the scintillation crystal

#averaging could be useful
def movingaverage(interval, window_size):
    window= np.ones(int(window_size))/float(window_size)
    return np.convolve(interval, window, 'same')
    

#Y=np.loadtxt('NIdata_5114_Co.txt') #for data stored in txt format
#Z=np.loadtxt('hxr_18730.txt') #for data from textronix osciloscope / old data
#Y=np.fromfile('Shotbin4',dtype='>f8', count=-1, sep='') #works for binaries
Y=np.fromfile('{}'.format(sh),dtype='>f8', count=-1, sep='') #works for binaries

#Z=movingaverage(Y,20)  #smoothing

kA=0.0394  # peak height to area coeficient
cshexp=np.exp(0.011*(U-789))   # scaling of peak height with power source voltage (exponential fit)
cshlin=np.exp(0.0125*(U-792))  # scaling of peak height with power source voltage (linearised fit)
print cshexp
Km=CsE/(cshexp*kA*mscint) #final coeficiet in keVs
KmJ=Km*1.6e-7 #final coefficient in nanoGreys
print KmJ
#kx=

X=np.linspace(0,len(Y)*(1/250000.0),len(Y)) #timeaxis generation
integ=scint.cumtrapz(Y,X) #integration of HXR signal

#plotting of HXR and and absorbed dose

fig, ax=plt.subplots(nrows=2)
ax[0].plot(X,Y, label='{} HXR'.format(sh))
ax[0].set_ylabel('Intensity [a.u.]')
ax[0].legend()
#print shape(X), shape(integ)
ax[1].plot(X[:-1],KmJ*integ, label='{} Dose'.format(sh), color='r')
ax[1].set_ylabel(r'$D$ [nGy]')
ax[1].set_xlabel(r'$t$ [ms]')
ax[1].legend()
print KmJ*integ[-1]

#plt.plot(np.linspace(0,len(integ)*(1/250000.0),len(integ)),k2*k*integ)
#plt.plot(X,Z)
#plt.title('HXR #18837, NI 250MHz, smoothed 10 S')
#plt.plot()
#plt.plot(K)
#plt.show()
#plt.savefig('Shot_HXR.png')

#saving and shoving the figure
plt.savefig('Golem_DOSE_{}.png'.format(sh)) 
plt.show()