Source code :: PlasmaPosition

[Return]
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
#!/usr/bin/env python
# -*- coding: utf-8 -*-


# ======================   GEFIT 0.1 (Golem EFIT) ======================================
# plama position calculation is based on the least square fit of the signal from mirnov coils a plasma current.
#A position and new plasma current are calculated from the furmula for loop wire is calculeted mag, filed in the coils and compared 
#with measured values. It should by robust and it was tested on the poor data from mirnov coils. 
#magnetic filed from external coils can be calculated by setting extCoilsPos and extCoilsCurr, however it was not tested yet. 
#
# Autor: Tomáš Odstrčil   tomasodstrcil@gmail.com

#TODO problémy: při velké disrupci dojde k tak masivnímu skoku že selže 
#integrace a nastane  integrační drift na některé z cívek. Udělá to chybu v signálu i 20%

import time
#t = time.time()
from numpy import *
from numpy.linalg import norm,lstsq
#from scipy.interpolate import interpolate
from scipy.optimize import minimize 
from multiprocessing import Process, Pool, cpu_count
#from Deconvolution import * 
from pygolem_lite.modules import deconvolveExp
import sys   
from MagFieldCalc import *
#print 'import time PP:', time.time()-t


#constants
a = 0.085		#[m]
R_0 = 0.4		#[m]



#http://www.phy.uct.ac.za/courses/python/examples/fitresonance.py
#http://www.scipy.org/doc/api_docs/SciPy.optimize.minpack.html

	
def directConductorApprox(B): 
    R =  (B[3]-B[1])/(B[3]+B[1])*-0.093
    Z  = (B[4]-B[2])/(B[4]+B[2])*-0.093
    r = hypot(Z,R) 
    if r > a:
	Z *= a/r
	R *= a/r
    R+= R_0
    return R,Z
    





def fun(x,param,B_coils,y,Bout):
    #print x.shape, param

    B = Bloop_analytic(x[:,0],x[:,1],param[0]+R_0,param[1],Bout)*param[2]
    B += B_coils    #NOTE musí být správná orientace proudů! meřených flukama
    B *=  x[:,2:]
   
    y[1:] = sqrt(sum(square(B),1))
    y[0] = abs(param[2])
    
    return y

	
def MainCycle((tvec,data,x,sigma,Bext,extCoilsCurr)):
    #print data.shape
    #data = double(data)
    #x = double(x)
    #tvec = double(tvec)
    #sigma = double(sigma)
    #Bext = double(Bext)
    #extCoilsCurr = double(extCoilsCurr)
    #print data.shape
    #exit()
    
    N = size(data,0)
    n_det = size(data,1)

    position = empty((N, 2))

    retrofit = empty((N,n_det))
    reziduum = empty(N)
    
    I0 = median(data[:,0])
    steps = 0   
    y = empty(size(x,0)+1)
    Bout = empty((size(x,0),2))
    for t in xrange(N): 
        #print t 
        #sys.stdout.write('\r'+str(t))
   
        #sys.stdout.flush()
        if all(isfinite(sigma)):
            R,Z = directConductorApprox(data[t,:])
            R+=0.04
        else:  R,Z =  R_0,0
	scale = 1/array((1,1,0.5/I0))   #1/a/2    
	params0 = array([R-R_0,Z,data[t,0]])

	if Bext == None:
	    B_coils = 0
	else:
	    B_coils = dot(Bext,extCoilsCurr[t,:])    
        
        T = time.time()
        
#def fitfun(params):
    
    #retrofit = fun(x,params*scale, B_coils,y,Bout)
    #retrofit-= data
    #retrofit/= sigma
    #retrofit *=retrofit
    #res = norm(retrofit)**2
    #res+= norm((params*scale)[:2]/(a*1.0))**50        
    ##res = norm((data[t,:] - fun(x,params*scale, B_coils,y,Bout))/sigma)**2+norm((params*scale)[:2]/(a*1.0))**50        
    
    #return res
        
	fitfun = lambda params : norm((data[t,:] - fun(x,params*scale, B_coils,y,Bout))/sigma)**2+norm((params*scale)[:2]/(a*1.0))**50	
       
	res = minimize(fitfun, params0/scale, method='Nelder-Mead',options={'xtol':1e-4, 'ftol' : 1e-2})
        #res = minimize(fitfun, params0*scale, method='BFGS',options={'xtol':1e-4, 'ftol' : 1e-2})

        #print '\n', time.time()-T
        #exit()
        #sys.stdout.write('\r'+str(time.time()-T))
   
        #sys.stdout.flush()
        
	if  not res.success:
	    res.x[:] = params0		
	    res.fun  = 1e6

	    print 'position calculation failure in ', str(tvec[t]), 's'
	else:
	    params0 =  res.x*scale

	
	position[t,:] = (res.x*scale)[:2]
	position[t,0]+= R_0
	steps+= res.nfev
	retrofit[t,:] = fun(x,params0,B_coils,y,Bout)
	reziduum[t] =  res.fun +sum(isinf(sigma))*1e4

    return (position,retrofit,reziduum)
    
#http://www.phy.uct.ac.za/courses/python/examples/fitresonance.py
#http://www.scipy.org/doc/api_docs/SciPy.optimize.minpack.html CalcPlasmaPosition(tvec_ds,detectorPos, detectorSignal_ds, detectorDriftError, Ip_ds,IpDriftError, shot_num = shot_num)

def CalcPlasmaPosition(tvec,detectorPos, detectorSignal,detectorDriftError,
                    Ip,IpDriftError,extCoilsPos = None,extCoilsCurr= None ):

#očekávám stejné časové rozdělení vstupního signálu, všechny položky jsou už absolutně kalibrované, jednotky SI
    print 'CalcPlasmaPosition'
    t1 = time.time()
    n =  len(tvec)

 
    coil_perpend=detectorPos-mean(detectorPos,axis=0)
    coil_perpend[:,1]*= -1   
    coil_perpend=coil_perpend[:,::-1]/hypot(coil_perpend[:,1],coil_perpend[:,0])[:,None]
   
    x = hstack((detectorPos,coil_perpend))

    data  = vstack((Ip, detectorSignal.T)).T
    sigma = hstack((IpDriftError, detectorDriftError))

    if extCoilsPos == None or extCoilsCurr == None:
	Bext = None
	
    else:
	Bext = zeros(shape(detectorPos)+(size(extCoilsPos,0),))

	for i in arange(size(extCoilsPos,0)):
	    Bext[...,i] = Bloop_analytic(x[:,0],x[:,1],extCoilsPos[i,0],extCoilsPos[i,0])

    
    #main calculation , multiprocess
    n_cpu = cpu_count()
    p = Pool(n_cpu)

    #split_ind =  array_split(range(n), n_cpu*4)
    #print data.shape
    #exit()
    data_spl = array_split(data, n_cpu*4,axis=0)
    tvec_spl = array_split(tvec, n_cpu*4)
    inputs = [(tvec_spl[i],data_spl[i],x,sigma,Bext,extCoilsCurr) for i in range(n_cpu*4)] 
    
    list_data = p.map(MainCycle,inputs )

    p.close()
    p.join()
    
    
    list_position = list()
    #list_position_dc = list()
    list_retrofit = list()
    list_reziduum = list()
    
    for (position,retrofit,reziduum) in list_data:
	list_position.append(position)
	#list_position_dc.append(position_dc)
	list_retrofit.append(retrofit)
	list_reziduum.append(reziduum)
	
    position = vstack(list_position)
    #position_dc = vstack(list_position_dc)
    retrofit = vstack(list_retrofit)
    reziduum = hstack(list_reziduum)


    radius = a-sqrt((position[:,0]-0.4)**2+position[:,1]**2)
    chi2 = mean((data-retrofit)**2,axis=0)/sigma**2

    print 'position calculated in %g s' % (time.time()-t1)
    print 'chi2: ', chi2,' total ', norm(chi2)
    
    return tvec, position, radius, reziduum, retrofit, data,chi2

  

    
def RemoveDriftsAuto(signal, Bt,Uloop,plasma_start,plasma_end,Bt_trigger,E_trigger, Bt_crosstalk=True, Uloop_crosstalk=False, Stabil_crosstalk=False, Trafo_effect=True ):
    n_det = size(signal, 1)-1

    index = 2
    t_max = min(Bt[-1,0],Uloop[-1,0],signal[-1,0])  # čas po který  to ty diagnostiky berou
    t_min = max(Bt[0,0],Uloop[0,0],signal[0,0])
    dt = mean(diff(signal[:,0]))
    signal = signal[:int(t_max/dt+1),:]
    tvec  = signal[:,0]
    n =  size(tvec)

    #resample to the same resolution
    if Uloop_crosstalk:
	i_Uloop = index
	index+=1
	Uloop = interp(tvec, Uloop[:,0],Uloop[:,1], left=0, right=None)

    if Bt_crosstalk:
	i_Bt = index
	index+=1
	Bt = interp(tvec, Bt[:,0],Bt[:,1], left=0, right=None)

    if Trafo_effect:
	t_exp = 21.2e-3 #[s]

	i_Traf_res = index
	Bt_conv,retrofit = deconvolveExp(Bt,-t_exp, dt,0,0)
	index+=1
    if Stabil_crosstalk:
	i_St = index
	index+=1

	Bt = interp(tvec, fluke[:,0],fluke[:,1], left=0, right=None)
	

    if isnan(plasma_start) or  isnan(plasma_end):	    
	plasma_start = Bt_trigger
	plasma_end = Bt_trigger
  
    #if removePlasmaDrift:
	#drift_error = mean(cumsum((diff(signal[:,1:], axis = 0))**2, axis = 0), axis = 1)
	#drift_error = hstack((0,drift_error))
	#plot(drift_error)
	#show()
	#plot(diff(signal[:,1:], axis = 0))
	#show()
	#i_drift = index
	#index+=1


    Bt_trigger_n = argmin(abs(tvec-Bt_trigger))

    
    E_trigger_n  = argmin(abs(tvec- E_trigger))
    plasma_start_n = max(argmin(abs(tvec-plasma_start)),E_trigger_n)
    plasma_end_n = max(argmin(abs(tvec-plasma_end)),E_trigger_n)
    

    interval0 = range(int(t_min/dt),Bt_trigger_n)   #clean area before any trigger
    interval1 = range(Bt_trigger_n+int(0.5e-3/dt), E_trigger_n)  # short time between Bt and CD trigger  
    interval2 = range(E_trigger_n,plasma_start_n)  #interval between  CD trigger   and breakdown
    interval3 = range(plasma_end_n,n )   # interval beween plasma end and torroidal field end tyristor
    interval_plasma = range(plasma_start_n,plasma_end_n)  #plasma 
    interval = interval1+interval2+interval3   #interval used for removing of the drift
    #print i
    mainBasis = zeros((n, index))
    if Bt_crosstalk:
	mainBasis[:,i_Bt] = Bt
    if Uloop_crosstalk:
	mainBasis[:,i_Uloop] = Uloop
    if Trafo_effect:
    	mainBasis[:,i_Traf_res] = Bt_conv
    if Stabil_crosstalk:
    	mainBasis[:,i_St] = fluke
    #if removePlasmaDrift:
	#mainBasis[:,i_drift] = drift_error

    #offset, integration drift
    offset = zeros(n)
    offset[interval0] = 1
    mainBasis[:,0] =  cumsum(offset, out = offset)
    offset = ones(n)
    offset[interval0] = 0
    mainBasis[:,1] =  cumsum(offset, out = offset)
    #offset = zeros(n)
    #offset[interval_plasma] = 1
    #mainBasis[:,2] =  cumsum(offset, out = offset)

    x,res,rank,s =  lstsq(mainBasis[interval,:],signal[interval,1:])

    corr_signal = signal[:,1:]-dot(mainBasis, x)
    res = sum((corr_signal[interval,:])**2,0)
    res*= 1e15/len(interval) 
    res = sqrt(res)
    
    #prepare graph
    frames = list()
    vlines = list()
    hlines = list()
    rectangles = list()
    vlines.append((plasma_end*1000,'--' ))
    vlines.append((plasma_start*1000,'--' ))
    hlines.append(( 0, '-.'))

    for i_det in arange(n_det):
	
	curves = list()
	curves.append((single(1e6*signal[:,i_det+1]),'Original signal', '--' ))
	curves.append((single(1e6*corr_signal[:,i_det]), 'Corrected signal', 'k-'))
	curves.append((single(1e6*dot(mainBasis[:,:2],x[:2,i_det])),'Integration offset','-' ))
	
	if Bt_crosstalk:
	    curves.append((single(1e6*mainBasis[:,i_Bt]*x[i_Bt,i_det]), 'Toroidal mag. field', '-'))	
	    
	if Trafo_effect:
	    curves.append((single(1e6*mainBasis[:,i_Traf_res]*x[i_Traf_res,i_det]),'Trafo mag. field','-'))
	
	if Uloop_crosstalk:
	    curves.append((single(1e6*mainBasis[:,i_Uloop]*x[i_Uloop,i_det]),'Uloop','-'))
	
	frames.append((1000*tvec,curves,'mc'+str(i_det*4+1)+' [ mV$\cdot$ms]','$\\chi^2$ = %2.1f'%res[i_det]))
	
	intervals = (interval1,interval2,interval3 )
    
    for inter in intervals:
	if len(inter) == 0:
	    continue
	(x_min,x_max) = (tvec[inter][0]*1e3, tvec[inter][-1]*1e3)
	rectangles.append((x_min,x_max))
    
    AutoRemoveGraph = (frames, vlines,hlines,rectangles)

    
    
    
    #save projection constants
    proj_file = 'mc1\tmc5\tmc9\tmc13\n'
    proj_file+='drift1 [V] '+str(x[0,:])+'\n'
    proj_file+='drift2 [V] '+str(x[1,:])+'\n'

    if Bt_crosstalk:
	    proj_file+='Bt crosstalk [V/T]\t'+str(x[i_Bt,:])+'\n'
    if Trafo_effect:
	    proj_file+='Trafo crosstalk (respons time %1.1f ms)[V/T]\t'%(1000*t_exp) +str(x[i_Traf_res,:])+'\n'
    if Uloop_crosstalk:
	    proj_file+='Uloop (residual) crosstalk [V/V]\t'+str(x[i_Uloop,:])+'\n'

    f = open('./constants/projections.txt', 'w')
    f.write(proj_file)
    f.close()
    
   

    return vstack((tvec,corr_signal.T)).T,res,AutoRemoveGraph
    
    
  

Navigation