##// END OF EJS Templates
schain
RSS
Update clase SpectralFitting
Alexander Valdez -
r1673:5478177ed526
parent child
Show More
Side by Side Unified
Context file:
r1673:5478177ed526 Loading diff...:
Show file before | Show file after | Hide comments
@@ -0,0 +1,132
1 import numpy
2
3 def setConstants(dataOut):
4 dictionary = {}
5 dictionary["M"] = dataOut.normFactor
6 dictionary["N"] = dataOut.nFFTPoints
7 dictionary["ippSeconds"] = dataOut.ippSeconds
8 dictionary["K"] = dataOut.nIncohInt
9
10 return dictionary
11
12 def initialValuesFunction(data_spc, constants):
13 #Constants
14 M = constants["M"]
15 N = constants["N"]
16 ippSeconds = constants["ippSeconds"]
17
18 S1 = data_spc[0,:]/(N*M)
19 S2 = data_spc[1,:]/(N*M)
20
21 Nl=min(S1)
22 A=sum(S1-Nl)/N
23 #x = dataOut.getVelRange() #below matches Madrigal data better
24 x=numpy.linspace(-(N/2)/(N*ippSeconds),(N/2-1)/(N*ippSeconds),N)*-(6.0/2)
25 v=sum(x*(S1-Nl))/sum(S1-Nl)
26 al1=numpy.sqrt(sum(x**2*(S1-Nl))/sum(S2-Nl)-v**2)
27 p0=[al1,A,A,v,min(S1),min(S2)]#first guess(width,amplitude,velocity,noise)
28 return p0
29
30 def modelFunction(p, constants):
31 ippSeconds = constants["ippSeconds"]
32 N = constants["N"]
33
34 fm_c = ACFtoSPC(p, constants)
35 fm = numpy.hstack((fm_c[0],fm_c[1]))
36 return fm
37
38 def errorFunction(p, constants, LT):
39
40 J=makeJacobian(p, constants)
41 J =numpy.dot(LT,J)
42 covm =numpy.linalg.inv(numpy.dot(J.T ,J))
43 #calculate error as the square root of the covariance matrix diagonal
44 #multiplying by 1.96 would give 95% confidence interval
45 err =numpy.sqrt(numpy.diag(covm))
46 return err
47
48 #-----------------------------------------------------------------------------------
49
50 def ACFw(alpha,A1,A2,vd,x,N,ippSeconds):
51 #creates weighted autocorrelation function based on the operational model
52 #x is n or N-n
53 k=2*numpy.pi/3.0
54 pdt=x*ippSeconds
55 #both correlated channels ACFs are created at the sametime
56 R1=A1*numpy.exp(-1j*k*vd*pdt)/numpy.sqrt(1+(alpha*k*pdt)**2)
57 R2=A2*numpy.exp(-1j*k*vd*pdt)/numpy.sqrt(1+(alpha*k*pdt)**2)
58 # T is the triangle weigthing function
59 T=1-abs(x)/N
60 Rp1=T*R1
61 Rp2=T*R2
62 return [Rp1,Rp2]
63
64 def ACFtoSPC(p, constants):
65 #calls the create ACF function and transforms the ACF to spectra
66 N = constants["N"]
67 ippSeconds = constants["ippSeconds"]
68
69 n=numpy.linspace(0,(N-1),N)
70 Nn=N-n
71 R = ACFw(p[0],p[1],p[2],p[3],n,N,ippSeconds)
72 RN = ACFw(p[0],p[1],p[2],p[3],Nn,N,ippSeconds)
73 Rf1=R[0]+numpy.conjugate(RN[0])
74 Rf2=R[1]+numpy.conjugate(RN[1])
75 sw1=numpy.fft.fft(Rf1,n=N)
76 sw2=numpy.fft.fft(Rf2,n=N)
77 #the fft needs to be shifted, noise added, and takes only the real part
78 sw0=numpy.real(numpy.fft.fftshift(sw1))+abs(p[4])
79 sw1=numpy.real(numpy.fft.fftshift(sw2))+abs(p[5])
80 return [sw0,sw1]
81
82 def makeJacobian(p, constants):
83 #create Jacobian matrix
84 N = constants["N"]
85 IPPt = constants["ippSeconds"]
86
87 n=numpy.linspace(0,(N-1),N)
88 Nn=N-n
89 k=2*numpy.pi/3.0
90 #created weighted ACF
91 R=ACFw(p[0],p[1],p[2],p[3],n,N,IPPt)
92 RN=ACFw(p[0],p[1],p[2],p[3],Nn,N,IPPt)
93 #take derivatives with respect to the fit parameters
94 Jalpha1=R[0]*-1*(k*n*IPPt)**2*p[0]/(1+(p[0]*k*n*IPPt)**2)+numpy.conjugate(RN[0]*-1*(k*Nn*IPPt)**2*p[0]/(1+(p[0]*k*Nn*IPPt)**2))
95 Jalpha2=R[1]*-1*(k*n*IPPt)**2*p[0]/(1+(p[0]*k*n*IPPt)**2)+numpy.conjugate(RN[1]*-1*(k*Nn*IPPt)**2*p[0]/(1+(p[0]*k*Nn*IPPt)**2))
96 JA1=R[0]/p[1]+numpy.conjugate(RN[0]/p[1])
97 JA2=R[1]/p[2]+numpy.conjugate(RN[1]/p[2])
98 Jvd1=R[0]*-1j*k*n*IPPt+numpy.conjugate(RN[0]*-1j*k*Nn*IPPt)
99 Jvd2=R[1]*-1j*k*n*IPPt+numpy.conjugate(RN[1]*-1j*k*Nn*IPPt)
100 #fft
101 sJalp1=numpy.fft.fft(Jalpha1,n=N)
102 sJalp2=numpy.fft.fft(Jalpha2,n=N)
103 sJA1=numpy.fft.fft(JA1,n=N)
104 sJA2=numpy.fft.fft(JA2,n=N)
105 sJvd1=numpy.fft.fft(Jvd1,n=N)
106 sJvd2=numpy.fft.fft(Jvd2,n=N)
107 sJalp1=numpy.real(numpy.fft.fftshift(sJalp1))
108 sJalp2=numpy.real(numpy.fft.fftshift(sJalp2))
109 sJA1=numpy.real(numpy.fft.fftshift(sJA1))
110 sJA2=numpy.real(numpy.fft.fftshift(sJA2))
111 sJvd1=numpy.real(numpy.fft.fftshift(sJvd1))
112 sJvd2=numpy.real(numpy.fft.fftshift(sJvd2))
113 sJnoise=numpy.ones(numpy.shape(sJvd1))
114 #combine arrays
115 za=numpy.zeros([N])
116 sJalp=zip(sJalp1,sJalp2)
117 sJA1=zip(sJA1,za)
118 sJA2=zip(za,sJA2)
119 sJvd=zip(sJvd1,sJvd2)
120 sJn1=zip(sJnoise, za)
121 sJn2=zip(za, sJnoise)
122 #reshape from 2D to 1D
123 sJalp=numpy.reshape(list(sJalp), [2*N])
124 sJA1=numpy.reshape(list(sJA1), [2*N])
125 sJA2=numpy.reshape(list(sJA2), [2*N])
126 sJvd=numpy.reshape(list(sJvd), [2*N])
127 sJn1=numpy.reshape(list(sJn1), [2*N])
128 sJn2=numpy.reshape(list(sJn2), [2*N])
129 #combine into matrix and transpose
130 J=numpy.array([sJalp,sJA1,sJA2,sJvd,sJn1,sJn2])
131 J=J.T
132 return J
Show file before | Show file after | Hide comments
@@ -0,0 +1,73
1 import numpy
2
3 def weightfit(w,year,doy,index,ht,beam):
4 if len(w) <= 0 : return w
5
6 if year == 2020 and doy == 265 :
7 if beam == 0 :
8 if index >= 120 and index <= 183 and ht >=30: w[0,100:120] =0
9
10 if beam == 1 :
11 if index >= 132 and index <= 180 and ht ==11:
12 w[0,90:120] = 0
13 if index >= 120 and index <= 183 and ht >=30:
14 w[0,100:120] = 0
15 if index >= 176 and index <= 176 and ht ==20: w[0,60:100] = 0.0
16 if index >= 177 and index <= 177 and ht >=24 and ht <=29: w[0,70:80] = 0.0
17 if index >= 177 and index <= 177 and ht >=30: w[0,70:80] = 0.0
18 if index >= 178 and index <= 178 and ht>=18: w[0,100:120] = 0.0
19 print('entra a la funcion w')
20 # If index eq 96 then stop ;stop 96
21
22 return w
23
24 def Vrfit(p0,year,doy,index,h,beam):
25 if year == 2020 and doy == 265 :
26 if beam == 0 :
27 if index == 177 and h==30 : p0[3]=20
28 if index == 184 and h>=27 : p0[3]=30
29 if index == 186 and h>=11 : p0[3]=30
30 if index == 190 and h==27 : p0[3]=27
31 if index == 191 and h==26 : p0[3]=28
32 if index == 191 and h>=35 : p0[3]=24
33 if index == 131 and h==29 : p0[3]=24
34 if index == 138 and h==11 : p0[3]=30
35 if index == 138 and h==14 : p0[3]=30
36 if index == 139 and h>=34 and h<=36 : p0[3]=20
37 if index == 142 and h>=34 : p0[3]=35
38 if index == 155 and h>=27 and h<=40: p0[3]=36
39 if index == 167 and h>=13 and h<=13: p0[3]=33
40 if index == 168 and h>=13 and h<=13: p0[3]=37
41 if index == 172 and h>=13 and h<=13: p0[3]=37
42 if index == 173 and h>=28 and h<=29: p0[3]=37
43 print('entra a la funcion vr')
44
45 if beam == 1 :
46 if index == 177 and h>=31 and h<=32: p0[3]=30
47 if index == 179 and h==27 : p0[3]=38
48 if index == 180 and h==23 : p0[3]=33
49 if index == 180 and h==26 : p0[3]=30
50 if index == 180 and h==29 : p0[3]=35
51 if index == 184 and h>=23 : p0[3]=30
52 if index == 185 and h==11 : p0[3]=20
53 if index == 185 and h>=29 : p0[3]=30
54 if index == 186 and h>=11 : p0[3]=30
55 if index == 190 and h==11 : p0[3]=28
56 if index == 191 and h>=33 : p0[3]=20
57 if index == 131 and h==11 : p0[3]=17
58 if index == 131 and h>=29 and h<=38 : p0[3]=18
59 if index == 138 and h>=31 : p0[3]=15
60 if index == 139 and h>=11 and h<=16 : p0[3]=25
61 if index == 140 and h==33 : p0[3]=26
62 if index == 142 and h==34 : p0[3]=20
63 if index == 149 and h>=28 and h<=32: p0[3]=28
64 if index == 163 and h>=11 and h<=11: p0[3]=30
65 if index == 167 and h>=11 and h<=16: p0[3]=33
66 if index == 167 and h>=34 and h<=40: p0[3]=36
67 if index == 169 and h>=26 and h<=49: p0[3]=33
68 if index == 172 and h>=28 and h<=28: p0[3]=37
69 if index == 173 and h>=20 and h<=38: p0[3]=37
70 if index == 174 and h>=25 and h<=25: p0[3]=34
71 if index == 175 and h>=27 and h<=35: p0[3]=30
72 #if h == 11: print('entra ',p0[3])
73 return p0
General Comments 0
You need to be logged in to leave comments. Login now