#!/usr/bin/python

# script takes as arguements: data_file pre std_num 
# script output:
# - prints transition temperature and standard deviation calculated using 'pre' as the number of data points to be used as the baseline, and std_num as the threshold number of standard deviations from the baseline considered to be phase separated
# - saves a plot of standard deviations from the baseline vs temperature for each ROI 

# data_file should be of the form:

#  "point	temp	StdDev(1)   StdDev(2)	StdDev(3)	StdDev(4)	StdDev(5)"	<--- headers are only for illustration, the input file should NOT contain any headers
#   1		60		26.633  	18.819  	30.897  	16.935  	30.993
#   2		60		26.636	    18.658  	30.772  	16.817  	31.13
#   3		60		26.641	    18.514  	30.897  	17.257      31.078
#   4		60		26.382	    18.749  	30.772  	16.838  	30.929
#   5		60		26.58	    18.762  	30.715  	16.82   	30.967
#   6		60		26.479  	18.696  	31.065  	16.978  	31.165
#   7		60		26.475	    18.431  	30.77   	16.74   	30.951
#   8		60		26.528  	18.582  	30.91   	16.758  	31.101
#   9		60		26.602  	18.638  	31.026  	16.767  	31.049
#   10		60		26.499  	18.883  	30.951  	16.72   	30.998
#   11		60		26.654  	18.848  	30.911  	16.565  	30.946
#   12		60		26.597  	18.688  	30.836  	16.758  	31.033
#   13		60		26.483  	18.706  	30.69   	16.647  	30.959
#   14		59.67	26.675  	18.792  	30.86   	16.811  	31.227
#   15		59.33	26.652  	18.842  	30.934  	16.79   	31.047
#   16		59.00	26.47   	18.726  	31.051  	16.857  	31.179
#   17		58.67	26.608  	18.882  	30.791  	16.89   	30.95
#   18		58.33	26.688  	18.771  	30.902  	16.751  	31.011
#   19		58.00	26.389  	18.762  	30.643  	16.777  	30.889
#   20		57.67	26.532  	18.768  	31.039  	16.761  	30.924
#   21		57.33	26.538  	18.683  	30.955  	16.806  	30.938
#   22		57.00	26.619  	18.859  	30.874  	16.811  	31.037
#   23		56.67	26.513  	18.658  	30.76   	16.604  	31.1
#   24		56.33	26.317  	18.88   	30.855  	16.756  	31.109
#   25		56.00	26.575  	18.739  	30.688  	16.712  	31.081
#   26		55.67	26.613  	18.685  	30.824  	16.668  	30.945
#   27		55.33	26.563  	18.777  	30.842  	16.731  	30.907
#   28		55.00	26.691  	18.734  	30.909  	16.864  	30.994
#   29		54.67	26.559  	18.855  	30.89	    16.536  	31.22
#   30		54.33	26.325  	18.841  	30.807  	16.643	    31.077
#   31		54.00	26.445  	18.747  	30.753  	16.641	    30.936

################################################################################################################################################################################################

import numpy as np
import matplotlib
import pylab as plt
import scipy
from sys import argv



script_name, data_file, pre, std_num = argv


data = np.loadtxt(data_file)	

temp = data[:,[1]]		

transition_temps = []	

pre = int(pre)					
std_num = int(std_num)			


def get_info(a):
	baseline_limit = np.mean(a[:pre]) + std_num * np.std(a[:pre])
	
	for x in a:
	    if x >= baseline_limit:
	       break		
	else:
	    x = None		

	array_loc = np.where(a == x)
	
	transition_temps.append(temp[array_loc])



get_info(data[:,[2]])		
get_info(data[:,[3]])
get_info(data[:,[4]])
get_info(data[:,[5]])
get_info(data[:,[6]])


Tp_avg = np.mean(transition_temps)
Tp_std = np.std(transition_temps)

print  np.around(Tp_avg, decimals = 2), np.around(Tp_std, decimals = 2)

def norm(a):						
	a = a - np.mean(a[pre])			

	return a						


col2 = norm(data[:,[2]])
col3 = norm(data[:,[3]])
col4 = norm(data[:,[4]])
col5 = norm(data[:,[5]])
col6 = norm(data[:,[6]])

means = sum([col2,col3,col4,col5,col6])/5

plt.rc('text', usetex=True)
plt.rc('font', family='serif', size=15)
plt.gcf().subplots_adjust(bottom=0.15)

plt.plot(temp, col2, 'k', alpha = 0.3)
plt.plot(temp, col3, 'k', alpha = 0.3)
plt.plot(temp, col4, 'k', alpha = 0.3)
plt.plot(temp, col5, 'k', alpha = 0.3)
plt.plot(temp, col6, 'k', alpha = 0.3)

plt.plot([50,0],[std_num, std_num],'r--', alpha = 0.4)


plt.minorticks_on()

plt.xlim(max(temp), min(temp)) # <- for LCST measurements, replace this line with this -> plt.xlim(min(temp), max(temp)) 

plt.ylim(min(col2), max(col2) * 1.2)

plt.xlabel(r'Temperature ($^\circ$C)')
plt.ylabel(r'$\Delta$ Pixel Intensity Standard Deviation')
plt.text(max(temp)-(max(temp)-min(temp))*0.1,max(col2)*1.1, "Transition Temperature = %.1f +/- %.1f$^\circ$C" %(Tp_avg,Tp_std), fontsize=15) # <- for LCST measurements, replace this line with this -> plt.text(min(temp)+(max(temp)-min(temp))*0.1,max(col2)*1.1, "Transition Temperature = %.1f +/- %.1f$^\circ$C" %(Tp_avg,Tp_std), fontsize=15)


plt.savefig(data_file[:-4] + '_plot.png', dpi = 1000)