modelSpectralFitting.py
132 lines
| 4.6 KiB
| text/x-python
|
PythonLexer
| r1396 | import numpy | |||
| def setConstants(dataOut): | ||||
| dictionary = {} | ||||
| dictionary["M"] = dataOut.normFactor | ||||
| dictionary["N"] = dataOut.nFFTPoints | ||||
| dictionary["ippSeconds"] = dataOut.ippSeconds | ||||
| dictionary["K"] = dataOut.nIncohInt | ||||
| return dictionary | ||||
| def initialValuesFunction(data_spc, constants): | ||||
| #Constants | ||||
| M = constants["M"] | ||||
| N = constants["N"] | ||||
| ippSeconds = constants["ippSeconds"] | ||||
| S1 = data_spc[0,:]/(N*M) | ||||
| S2 = data_spc[1,:]/(N*M) | ||||
| Nl=min(S1) | ||||
| A=sum(S1-Nl)/N | ||||
| #x = dataOut.getVelRange() #below matches Madrigal data better | ||||
| x=numpy.linspace(-(N/2)/(N*ippSeconds),(N/2-1)/(N*ippSeconds),N)*-(6.0/2) | ||||
| v=sum(x*(S1-Nl))/sum(S1-Nl) | ||||
| al1=numpy.sqrt(sum(x**2*(S1-Nl))/sum(S2-Nl)-v**2) | ||||
| p0=[al1,A,A,v,min(S1),min(S2)]#first guess(width,amplitude,velocity,noise) | ||||
| return p0 | ||||
| def modelFunction(p, constants): | ||||
| ippSeconds = constants["ippSeconds"] | ||||
| N = constants["N"] | ||||
| fm_c = ACFtoSPC(p, constants) | ||||
| fm = numpy.hstack((fm_c[0],fm_c[1])) | ||||
| return fm | ||||
| def errorFunction(p, constants, LT): | ||||
| J=makeJacobian(p, constants) | ||||
| J =numpy.dot(LT,J) | ||||
| covm =numpy.linalg.inv(numpy.dot(J.T ,J)) | ||||
| #calculate error as the square root of the covariance matrix diagonal | ||||
| #multiplying by 1.96 would give 95% confidence interval | ||||
| err =numpy.sqrt(numpy.diag(covm)) | ||||
| return err | ||||
| #----------------------------------------------------------------------------------- | ||||
| def ACFw(alpha,A1,A2,vd,x,N,ippSeconds): | ||||
| #creates weighted autocorrelation function based on the operational model | ||||
| #x is n or N-n | ||||
| k=2*numpy.pi/3.0 | ||||
| pdt=x*ippSeconds | ||||
| #both correlated channels ACFs are created at the sametime | ||||
| R1=A1*numpy.exp(-1j*k*vd*pdt)/numpy.sqrt(1+(alpha*k*pdt)**2) | ||||
| R2=A2*numpy.exp(-1j*k*vd*pdt)/numpy.sqrt(1+(alpha*k*pdt)**2) | ||||
| # T is the triangle weigthing function | ||||
| T=1-abs(x)/N | ||||
| Rp1=T*R1 | ||||
| Rp2=T*R2 | ||||
| return [Rp1,Rp2] | ||||
| def ACFtoSPC(p, constants): | ||||
| #calls the create ACF function and transforms the ACF to spectra | ||||
| N = constants["N"] | ||||
| ippSeconds = constants["ippSeconds"] | ||||
| n=numpy.linspace(0,(N-1),N) | ||||
| Nn=N-n | ||||
| R = ACFw(p[0],p[1],p[2],p[3],n,N,ippSeconds) | ||||
| RN = ACFw(p[0],p[1],p[2],p[3],Nn,N,ippSeconds) | ||||
| Rf1=R[0]+numpy.conjugate(RN[0]) | ||||
| Rf2=R[1]+numpy.conjugate(RN[1]) | ||||
| sw1=numpy.fft.fft(Rf1,n=N) | ||||
| sw2=numpy.fft.fft(Rf2,n=N) | ||||
| #the fft needs to be shifted, noise added, and takes only the real part | ||||
| sw0=numpy.real(numpy.fft.fftshift(sw1))+abs(p[4]) | ||||
| sw1=numpy.real(numpy.fft.fftshift(sw2))+abs(p[5]) | ||||
| return [sw0,sw1] | ||||
| def makeJacobian(p, constants): | ||||
| #create Jacobian matrix | ||||
| N = constants["N"] | ||||
| IPPt = constants["ippSeconds"] | ||||
| n=numpy.linspace(0,(N-1),N) | ||||
| Nn=N-n | ||||
| k=2*numpy.pi/3.0 | ||||
| #created weighted ACF | ||||
| R=ACFw(p[0],p[1],p[2],p[3],n,N,IPPt) | ||||
| RN=ACFw(p[0],p[1],p[2],p[3],Nn,N,IPPt) | ||||
| #take derivatives with respect to the fit parameters | ||||
| 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)) | ||||
| 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)) | ||||
| JA1=R[0]/p[1]+numpy.conjugate(RN[0]/p[1]) | ||||
| JA2=R[1]/p[2]+numpy.conjugate(RN[1]/p[2]) | ||||
| Jvd1=R[0]*-1j*k*n*IPPt+numpy.conjugate(RN[0]*-1j*k*Nn*IPPt) | ||||
| Jvd2=R[1]*-1j*k*n*IPPt+numpy.conjugate(RN[1]*-1j*k*Nn*IPPt) | ||||
| #fft | ||||
| sJalp1=numpy.fft.fft(Jalpha1,n=N) | ||||
| sJalp2=numpy.fft.fft(Jalpha2,n=N) | ||||
| sJA1=numpy.fft.fft(JA1,n=N) | ||||
| sJA2=numpy.fft.fft(JA2,n=N) | ||||
| sJvd1=numpy.fft.fft(Jvd1,n=N) | ||||
| sJvd2=numpy.fft.fft(Jvd2,n=N) | ||||
| sJalp1=numpy.real(numpy.fft.fftshift(sJalp1)) | ||||
| sJalp2=numpy.real(numpy.fft.fftshift(sJalp2)) | ||||
| sJA1=numpy.real(numpy.fft.fftshift(sJA1)) | ||||
| sJA2=numpy.real(numpy.fft.fftshift(sJA2)) | ||||
| sJvd1=numpy.real(numpy.fft.fftshift(sJvd1)) | ||||
| sJvd2=numpy.real(numpy.fft.fftshift(sJvd2)) | ||||
| sJnoise=numpy.ones(numpy.shape(sJvd1)) | ||||
| #combine arrays | ||||
| za=numpy.zeros([N]) | ||||
| sJalp=zip(sJalp1,sJalp2) | ||||
| sJA1=zip(sJA1,za) | ||||
| sJA2=zip(za,sJA2) | ||||
| sJvd=zip(sJvd1,sJvd2) | ||||
| sJn1=zip(sJnoise, za) | ||||
| sJn2=zip(za, sJnoise) | ||||
| #reshape from 2D to 1D | ||||
| sJalp=numpy.reshape(list(sJalp), [2*N]) | ||||
| sJA1=numpy.reshape(list(sJA1), [2*N]) | ||||
| sJA2=numpy.reshape(list(sJA2), [2*N]) | ||||
| sJvd=numpy.reshape(list(sJvd), [2*N]) | ||||
| sJn1=numpy.reshape(list(sJn1), [2*N]) | ||||
| sJn2=numpy.reshape(list(sJn2), [2*N]) | ||||
| #combine into matrix and transpose | ||||
| J=numpy.array([sJalp,sJA1,sJA2,sJvd,sJn1,sJn2]) | ||||
| J=J.T | ||||
| return J | ||||
