#!/usr/bin/env python # -*- coding: utf-8 -*- ########################################################################################################### ############# Algorithm which from ion intensities calculate the coefitients to estimate electron t ############# temperature and impurity density. ########################################################################################################### ion_names = ['HI','OI','OII', 'OIII', 'HeI', 'CII', 'CIII', 'NI','NII', 'NIII','M1', 'MIV' ] element_names = ['H','O','He','C','N','M1','M2'] ions_index = [[0], [1,2,3], [4], [5,6], [7,8,9], [10],[11]] single_ion = [0,2,5,6] multi_ion = [1,3,4] import matplotlib #matplotlib.rcParams['backend'] = 'Agg' #matplotlib.rc('font', size='10') #matplotlib.rc('text', usetex=True) # FIXME !! nicer but slower !!! from numpy import * from matplotlib.pyplot import * import os from numpy.linalg import norm from scipy.linalg import qr,inv,pinv,eig,cholesky,solve_triangular from scipy.stats.stats import pearsonr,hmean from scipy.sparse import block_diag from scipy.optimize import fmin_powell from numpy.matlib import repmat from time import time from scipy.stats.stats import nanmean,nanmedian n_ions = 12 n_elemets = 7 k = 2 #polynome order err = load('err.npy') T = load('temperature.npy') shots = load('shots.npy') n = len(T) errT = err[-(1+n_elemets)*n:-n_elemets*n] #errDensity = reshape(err[-n*n_elemets:],(n,n_elemets)).T T_spitzer = zeros_like(T) T_spitzer_err = zeros_like(T) T_all= zeros(11000-7100) shots_all = arange(7100,11000) for shot in range(7100,11000): #print len(where(shot == int_(shots)) shot_index = where(abs(floor(shots)- shot) < 0.5)[0] n_spec = size(shot_index) ##if n_spec == 0 : ##continue try: #print 'loadim', shot try: data = load('./electronTemperatures/'+str(shot)+'/electron_temperature.npz') temp = data['data']*data['scale'] #plot(temp) #show() tvec= linspace(data['t_start'], data['t_end'], len(temp)) try: plasma_start = loadtxt('./electronTemperatures/'+str(shot)+'/PlasmaStart') plasma_end = loadtxt('./electronTemperatures/'+str(shot)+'/PlasmaEnd') except: plasma_start = amin(tvec[isfinite(temp)]) plasma_end = amax(tvec[isfinite(temp)]) #print plasma_start,plasma_end #print 'nacteno npz', (temp) except: #ss temp = loadtxt('./electronTemperatures/'+str(shot)+'/ElectronTemperatureMedianFilter.txt') tvec = temp[:,0] temp = temp[:,1] savez('./electronTemperatures/'+str(shot)+'/electron_temperature.npz', data=temp, scale=1, t_start=tvec[0],t_end=tvec[-1]) #print shape(temp) #try: plasma_start = loadtxt('./electronTemperatures/'+str(shot)+'/PlasmaStart') plasma_end = loadtxt('./electronTemperatures/'+str(shot)+'/PlasmaEnd') #except: #print 'no lasma start' #continue except: print shot, ' failured' continue print shot if shot > 9300: temp*= 2 else: tvec/= 1000 plasma_start-= 0.5 plasma_start/= 1000 plasma_end/=1000 #plot(tvec,temp) #xlim(plasma_start,plasma_end) #ylim(0,80) #savefig('./electronTemperatures/'+str(shot)+'.png') #clf() T_all[shot-7100] = median(temp[(tvec > plasma_start) * (tvec< plasma_end)]) tvec_spec = linspace(plasma_start,plasma_end,n_spec+1) #plot(temp[:,0],temp[:,1]) #xlim(plasma_start,plasma_end) #ylim(0,100) #show() for i in range(n_spec): shot_temp = temp[(tvec > tvec_spec[i]) * (tvec< tvec_spec[i+1])] T_spitzer[shot_index[i]] = nanmedian(shot_temp) T_spitzer_err[shot_index[i]] = nanmean(abs(shot_temp-nanmedian(shot_temp)))*1.4 #print T_spitzer[shot_index] #exit() plot(shots_all,T_all, '.') show() #plot(T_spitzer,-T, '.',markersize = 0.5) ##errorbar(T_spitzer,-T,yerr = errT, xerr = T_spitzer_err, fmt='.',capsize = 0,linewidth = 0.1,markersize = 0.5) #ylim(-4,2) #xlim(3,100) #show() #plot(shots,T_spitzer) ind = where((T_spitzer!= 0)) T = -T[ind] errT = errT[ind] T_spitzer = T_spitzer[ind] T_spitzer_err = T_spitzer_err[ind] n = size(ind) print n iter = 0 def f(var,x,xerr,y,yerr): #print sum(x),norm((xerr)) T_spitz = var[-n:] a = var[0] b = var[1] c = var[2] d = var[3] #y = param[0] #yerr = param[1] #x = param[2] #xerr = param[3] #T_new = zeros_like(T_spitz) #print shape(T_new), shape(y), shape(yerr), shape(T_spitz), shape(x), shape(xerr) #ind = T_spitz*a+b> 1e-3 #T_new[ind] = log((T_spitz*a+b)[ind]) T_new = d*arctan(T_spitz*a+b)+c #T_new = boundary = 0 #ind = argsort(T_spitz) #print shape(ind) #boundary = norm(diff(T_new[ind]) < 0) #boundary = 0#norm((T_spitz*a+b)[~ind]) #print shape(T_new), shape(y), shape(yerr), shape(T_spitz), shape(x), shape(xerr) #resid = norm((T_new-y)/yerr)**2+norm((T_spitz-x)/xerr)**2 Resid = hstack(((T_new-y)/yerr, (T_spitz-x)/xerr)) max_Resid = double(amax(abs(Resid))) Resid/= max_Resid Resid **= 2 chi2 = max_Resid*sum(Resid/sqrt((1/max_Resid)**2 + Resid))/(size(y)+size(x)-size(var)) #print norm((T_new-y)/yerr),norm((x)), sum(x),norm((1/xerr)) #chi2 = resid/(size(y)+size(x)-size(var)) #print resid global iter if iter%1000 ==0: print iter, chi2,boundary iter+=1 #return return chi2+boundary def TotalLqr(f, x,xerr,y,yerr): #plot(yerr) #show() #plot(x) #show() xerr[~isfinite(xerr)] = infty #print where(1-isfinite(xerr)) #xerr+= 0.01 #print sum(x)0.231930954473 -2.25248138786 -0.934499469238 x0 = hstack((0.23,-2.25,-0.93,2,x)) var = x0 try: #aa var = load('var.npy') except: var = fmin_powell(f,x0,args=(x,xerr,y,yerr), xtol=1e-5,ftol=1e-5,maxfun=1e6,maxiter=1e7) save('var',var) T_spitz = var[-n:] a = var[0] b = var[1] c = var[2] d = var[3] #yerr += 0.1 #xerr += 1 print a, b ,c,d errorbar(x,y,yerr=yerr, xerr=xerr, fmt='.',capsize = 0,linewidth = 0.1,markersize = 0.5) T_new = d*arctan(T_spitz*a+b)+c plot(T_spitz,T_new, '.') show() plot(x,y,'.',markersize = 0.5) T_new = d*arctan(T_spitz*a+b)+c plot(T_spitz,T_new, '.') show() ind = where((y-c)/d > pi) y[ind] = a*0.95*pi+c T = (tan((y-c)/d)-b)/a Terr = 1/(a*d)/cos((y-c)/d)**2*yerr errorbar(x,T,yerr=Terr,xerr=xerr,fmt='.',capsize = 0,linewidth = 0.1,markersize = 0.5) show() plot(x,T,'.',markersize = 0.5) plot([0,0],[60,60],'-') xlabel('T spectrometer [eV]') ylabel('T spitzer [eV]') ylim(0,60) xlim(0,60) show() errorbar(arange(n),T,yerr=Terr,capsize = 0,linewidth = 0.1,markersize = 0.5) show() print '-------------------' TotalLqr(f, T_spitzer,T_spitzer_err,T,errT) #1.21772166429