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. #-----------------------------
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.

Recent Pastes

Raw Paste

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