参考文献
铬稳定同位素双稀释剂法分析技术与应用_李理
没有深究原理,只是写个代码.猜测其他双稀释剂方法也可以用类似代码.
Python写的,方便使用.
from math import log,pow #define variables while(True): Debug = False f1 = 0.8 f2 = 1.2 M50 = 49.946046 M52 = 51.940509 M53 = 52.940651 M54 = 53.938882 Rs50 = 21.0956 Rs53 = 0.087417 Rs54 = 21.5006 Rn50 = [0 for i in range(100)] Rn53 = [0 for i in range(100)] Rn54 = [0 for i in range(100)] # Rn50[0] = (0.05186) # Rn53[0] = (0.11339) # Rn54[0] = (0.02822) Rn50[0] = (0.05186) Rn53[0] = (0.11339) Rn54[0] = (0.02822) Rm50 = [0 for i in range(12)] Rm53 = [0 for i in range(12)] Rm54 = [0 for i in range(12)] Rs = [0 for i in range(12)] Pn = [0 for i in range(100)] p = 0 N = 100 print("------------------------------------------------------------------------------") print("Double Spike Calibration Program For chromium(Cr) Isotopes -------------------") print("Author:Kai Li Contact:[email protected]") print("Only For Reasearch Use--------------------------------------------------------") print("Spike Isotope Value is,Rs50 = 21.0956 | Rs53 = 0.087417 | Rs54 = 21.5006 -----") print("------------------------------------------------------------------------------") if Debug: Rm50[0] = (2.201705) Rm53[0] = (0.110168) Rm54[0] = (2.197483) N = 100 else: print("Measured Data(eg:Rm50=2.201705 Rm53=0.110168 Rm54=2.197483 )") Rm50[0] =(float(input("Please Enter Rm50:"))) Rm53[0] =(float(input("Please Enter Rm53:"))) Rm54[0] =(float(input("Please Enter Rm54:"))) #N = (int(input("Please Enter Iter times:"))) for j in range(N): for i in range(11): Rs[i] = ( Rs50 * (Rn50[j] -Rm50[i]) * (Rm54[i] - Rs54)) / ( Rs54 * (Rn54[j] - Rm54[i]) * (Rm50[i] - Rs50)) p = (log( Rs[i] / (Rs50 / Rs54)) / log( M50/M54) ) * f1 Rm50 [i+1] = Rm50[i] / pow( M50/M52, p) Rm53 [i+1] = Rm53[i] / pow( M53/M52, p) Rm54 [i+1] = Rm54[i] / pow( M54/M52, p) Rm50[0] = Rm50[11] Rm53[0] = Rm53[11] Rm54[0] = Rm54[11] Rn53[j+1] = ( Rm53[0] - Rs53)*(Rm50[0] - Rn50[j])/(Rs50 - Rm50[0]) + Rm53[0] if j==0: Pn[j] = ( log( Rn53[j+1] / Rn53[j]) / log( M53/M52) ) * f2 elif j==N-1: Pn[j] = ( log( Rn53[j+1] / Rn53[j]) / log( M53/M52) ) else: Pn[j] = ( log( Rn53[j+1] / Rn53[j]) / log( M53/M52) ) * f2 - (f2-1)*Pn[j-1] Rn50[j+1] = Rn50[j] * pow(M50/M52, Pn[j]) Rn53[j+1] = Rn53[j] * pow(M53/M52, Pn[j]) Rn54[j+1] = Rn54[j] * pow(M54/M52, Pn[j]) if (Rn50[j+1] - Rn50[j])<0.0000001 and (Rn53[j+1] - Rn53[j])<0.0000001 and (Rn54[j+1] - Rn54[j])<0.0000001: break elif (j+1) == N: print("input N again") print(j, end='') print(":","Cr50/Cr52=",round(Rn50[j+1],7), "|Cr53/Cr52=", round(Rn53[j+1],7), "|Cr54/Cr52=", round(Rn54[j+1],7)) print("") input("Press to Continue")