利用氢氘比和有效z确定密度剖面

defaultVars(Zeff=2,fhd=0.05)
# IONS_1: D; IONS_2: H;  IONS_3: impurity
# 94249 Zeff=2.1,fhd=0.111111
# 93320 Zeff=1.2,fhd=0.05
# 94238 Zeff=2.5,fhd=0.3/0.7
Zeff=4
fhd=0.1111
profiles=OMFIT['new_kineticEFIT']['profiles']['CASE']['94249']['PROFILES']
profiles['fHD']=fhd
ne=profiles['ne']*1.
ni=profiles['ni']
profiles['Zeff']=profiles['Zeff']*0+Zeff

qH=profiles['IONS_2']['charge']
qD=profiles['IONS_1']['charge']
qC=profiles['IONS_3']['charge']

Zc=profiles['IONS_3']['mass']
ZH=profiles['IONS_2']['mass']
ZD=profiles['IONS_1']['mass']

A = fhd*qH+qD

B = ZH*Zeff*fhd+ZD*Zeff-ZH**2*fhd-ZD**2

C = Zc**2-Zc*Zeff

#D = A*Zc**2-B*qC

profiles['IONS_3']['density']= ne/A/(C/B+qC/A)

profiles['IONS_1']['density']= C/B*profiles['IONS_3']['density']

profiles['IONS_2']['density']= profiles['IONS_1']['density']*fhd

profiles['ni']=profiles['IONS_1']['density']

for i in range(3):
    node=profiles['IONS_'+str(i+1)]
    node['ni_over_ne']=node['density']/ne

figure()
plot(ne,'-',label='ne')
plot(profiles['IONS_1']['density'],label='ni1')
plot(profiles['IONS_2']['density']/fhd,label='ni2/fhd')
plot(profiles['IONS_3']['density'],label='ni3')
plot(profiles['ni'],'--',label='ni')
legend(loc=1)