PYTHON 23
DatafitSOL.py Guest on 2nd May 2021 02:27:28 AM
  1. #-------------------------
  2. # datafitSOL.py        
  3. # read SOL data file
  4. #------------------------
  5. import matplotlib.pyplot as plt
  6. from numpy import *
  7. #------------------------------------------
  8. def linfit(x,y,wt):
  9. # Function to fit a line
  10. # Inputs -
  11. #   x - Independent variable
  12. #   y - Dependent variable
  13. #   wt - weight
  14. # Outputs -
  15. #   a_fit - Fit parameters; a(1) is intercept, a(2) is slope
  16. #   sig_a - Estimated error in the parameters a()
  17. #-----------------------------------------
  18.     N =   size(x)
  19.     s =   sum(wt)              
  20.     sx=dot(x,wt)
  21.     sy=dot(y,wt)
  22.     sxy=0.0;  sxx=0.0
  23.     for k in num:
  24.         sxy=sxy+x[k]*y[k]*wt[k]
  25.         sxx=sxx+x[k]*x[k]*wt[k]
  26. # % matrix method
  27.     alpha = [[s, sx], [ sx, sxx]];
  28.     betas  = [sy, sxy];
  29. #
  30.     errmat = linalg.inv(alpha);
  31.     param  = dot(betas,errmat);
  32. #  the fitted parameters
  33.     a_fit = param
  34. #    print ' a_fit  ' %a_fit
  35. #    return a_fit
  36. #  calculate errors from error matrix
  37.     sig_a0 = sqrt(errmat[[0],[0]])
  38.     sig_a1 = sqrt(errmat[[1],[1]])
  39.     sig_a=[0,0]
  40.     sig_a[0]=sig_a0
  41.     sig_a[1]=sig_a1
  42. #    print 'sig ', sig_a[0],sig_a[1]
  43.     return a_fit,sig_a
  44. #---------------------------------------
  45. #   calculate chi sq ----------------
  46. def chi2(y,yy,wtt):
  47.     ylen=len(y)
  48.     chs=0.0
  49.     for iy in range(0,ylen):
  50.         chs=chs+((y[iy]-yy[iy])**2)*wtt[iy]
  51. #    print 'chi sq %8.2f' %chs
  52.     return chs
  53. #-----------------------------
  54. # read data file
  55. num=[] ; x=[]; errdat=[]; errors=[]; wtt=[]
  56. #f=open("lindata.txt")
  57. f=open("linsol.txt")
  58. k=0
  59. for line in f:
  60.     win=line.split()
  61. #    print win
  62.     xin=float(win[0])
  63.     datin=float(win[1])
  64.     erin=float(win[2])
  65.     x.append(xin)
  66.     errdat.append(datin)
  67.     errors.append(erin)
  68. # weight factor
  69.     wtt.append((1.0/erin)**2)
  70.     num.append(k)
  71.     k=k+1
  72. #--------------------------------
  73. #for k in num:
  74. #    print '%d  %4.2f  %6.2f  %6.2f  %6.4f' %(k,x[k],errdat[k],errors[k],wtt[k])
  75. #----------------------------------
  76. # fit the line
  77. a_fit,sig_a=linfit(x,errdat,wtt)
  78. print "a_fit %6.3f  %10.4e" %(a_fit[0],a_fit[1])
  79. print "sig_a %6.3f  %10.4e" %(sig_a[0],sig_a[1])
  80. #-----------------------------------
  81. # calculate chi sq for best fit
  82. xlen=len(x)
  83. yy=zeros(xlen)
  84. inter=a_fit[0]
  85. slpe=a_fit[1]
  86. for iiy in range(0,xlen):
  87.     yy[iiy]=inter + slpe*x[iiy]
  88. chs=chi2(errdat,yy,wtt)
  89. print ' %5d points      chi sq %6.2f' %(xlen,chs)
  90. #-----------------------------------
  91. # chi sq as a function of slope
  92. slope=zeros(80); chst=zeros(80)
  93. yyy=zeros(xlen)
  94. for k in range(0,79):
  95.     slope[k]=0.00015+0.0000005*k
  96.     for iiy in range(0,xlen):
  97.         yyy[iiy]=inter+ slope[k]*x[iiy]
  98.     chst[k]=chi2(errdat,yyy,wtt)
  99. #-----------------------------------------
  100. plt.figure()
  101. # ---------  plot chi sq as a function of slope
  102. #plt.subplot(211)
  103. plt.plot(slope,chst,lw=0)
  104. plt.errorbar(slope,chst,0,fmt='o')
  105. plt.xlabel("slope")
  106. plt.ylabel("chisq")
  107. plt.ylim(chs-1,chs+4)
  108. plt.xlim(0.00018,0.00019)
  109. #plt.xlim(slpe-4,slpe+4)
  110. plt.grid(True)
  111. #plt.title("%4d points   chisq %6.2f    slope %7.2f     inter %7.2f " %(xlen,chs,a_fit[1],a_fit[0]))
  112. #-------- plot data and fitted line -----------
  113. plt.figure()
  114. #plt.subplot(212)
  115. plt.plot(x,errdat,lw=0)
  116. plt.plot(x,yy,lw=1)
  117. plt.errorbar(x,errdat,errors,fmt='o')
  118. plt.xlabel("x")
  119. plt.ylabel("data")
  120. plt.ylim(9.5,10.5)
  121. #plt.ylim(0,1100)
  122. plt.grid(True)
  123. plt.title("%4d points    chisq %6.2f    slope %9.3e   intercept %6.3f " %(xlen,chs,a_fit[1],a_fit[0]))
  124. plt.savefig("SOLdata.png")
  125. plt.show()

Paste-bin is for source code and general debugging text.

Login or Register to edit, delete and keep track of your pastes and more.

Raw Paste

Login or Register to edit or fork this paste. It's free.