铬同位素双稀释剂矫正代码python(geokai)

参考文献

铬稳定同位素双稀释剂法分析技术与应用_李理

没有深究原理,只是写个代码.猜测其他双稀释剂方法也可以用类似代码.

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:geokai@126.com----------------------------------------")
    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")