##// END OF EJS Templates
schain
RSS
Updates for EW-Drifts
Percy Condor -
r1382:e1808c191132
parent child
Show More
Side by Side Unified
Context file:
r1382:e1808c191132 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
General Comments 0
You need to be logged in to leave comments. Login now