1 | from model import Model |
---|
2 | from pySAXS.LS.LSsca import Qlogspace |
---|
3 | from pySAXS.LS.LSsca import * |
---|
4 | import numpy |
---|
5 | |
---|
6 | class Pcylmultimin(Model): |
---|
7 | ''' |
---|
8 | Multilayer Cylinder |
---|
9 | by OT 10/06/2009 |
---|
10 | ''' |
---|
11 | def PcylmultiminFunction(self,q,par): |
---|
12 | R=[par[0],par[1],par[2],par[3],par[4],par[5]] |
---|
13 | rho=[par[6],par[7],par[8],par[9],par[10],par[11],par[12]] |
---|
14 | L=par[13] |
---|
15 | rhos=par[14] |
---|
16 | return self.Pcylmulti(q,R,rho,L,rhos) |
---|
17 | |
---|
18 | def Pcylmulti(self,q,R,rho,L,rhos): |
---|
19 | """ |
---|
20 | This function calculates the P(q) of a cylinder of length L and and multilayers R of density Rho the last one being solvent |
---|
21 | """ |
---|
22 | Pcylnew=1.*q |
---|
23 | Pcylold=1.*q |
---|
24 | for i in range(len(q)): |
---|
25 | N=1. |
---|
26 | alpha0=1e-12 |
---|
27 | alphaN=pi/2 |
---|
28 | |
---|
29 | #---------------------------------------------------------------------------------------------------- |
---|
30 | """Somme sur les ecorces""" |
---|
31 | dPeco=0. |
---|
32 | for j in range(len(R)): |
---|
33 | dPeco=dPeco+(rho[j]-rho[j+1])*(2.*R[j]*SPspecial.j1(q[i]*R[j]*numpy.sin(alpha0))/(q[i]*numpy.sin(alpha0))) |
---|
34 | dPcyl0=numpy.sin(alpha0)*(numpy.sin(q[i]*L*numpy.cos(alpha0)/2.)/(q[i]*L*numpy.cos(alpha0)/2.)*dPeco)**2 |
---|
35 | |
---|
36 | #---------------------------------------------------------------------------------------------------- |
---|
37 | dPeco=0. |
---|
38 | for k in range(len(R)): |
---|
39 | dPeco=dPeco+(rho[k]-rho[k+1])*(2.*R[k]*SPspecial.j1(q[i]*R[k]*numpy.sin(alphaN))/(q[i]*numpy.sin(alphaN))) |
---|
40 | dPcylN=numpy.sin(alphaN)*(numpy.sin(q[i]*L*numpy.cos(alphaN)/2.)/(q[i]*L*numpy.cos(alphaN)/2.)*dPeco)**2 |
---|
41 | #---------------------------------------------------------------------------------------------------- |
---|
42 | |
---|
43 | dPcylold=(dPcyl0+dPcylN)/2. |
---|
44 | Pcylold[i]=pi/2*dPcylold |
---|
45 | test=1. |
---|
46 | while(abs(test)>1e-4): |
---|
47 | N=2.*N |
---|
48 | dalpha=(pi/2.-alpha0)/N |
---|
49 | alpha=numpy.arange(dalpha,pi/2.,2.*dalpha) |
---|
50 | #---------------------------------------------------------------------------------------------------- |
---|
51 | dPeco=0.*alpha |
---|
52 | for l in range(len(R)): |
---|
53 | dPeco=dPeco+(rho[l]-rho[l+1])*(2.*R[l]*SPspecial.j1(q[i]*R[l]*numpy.sin(alpha))/(q[i]*numpy.sin(alpha))) |
---|
54 | dPcyl=numpy.sin(alpha)*(numpy.sin(q[i]*L*numpy.cos(alpha)/2.)/(q[i]*L*numpy.cos(alpha)/2.)*dPeco)**2 |
---|
55 | #---------------------------------------------------------------------------------------------------- |
---|
56 | dPcylnew=sum(dPcyl)+dPcylold |
---|
57 | dPcylold=dPcylnew |
---|
58 | Pcylnew[i]=dalpha*dPcylnew |
---|
59 | test=(Pcylnew[i]-Pcylold[i])/Pcylold[i] |
---|
60 | Pcylold[i]=Pcylnew[i] |
---|
61 | return Pcylold*numpy.pi**2*L**2*1e-48*rhos |
---|
62 | |
---|
63 | ''' |
---|
64 | parameters definition |
---|
65 | Model(10,Pcylmultimin,Qlogspace(1e-4,1.,500.), |
---|
66 | ([13.,14.,15.,16.,17.,18.,9.418e+010,8.798e+010,4.041e+011,3.925e+011,4.224e+011,3.902e+011,9.418e+010,400.,1e17]), |
---|
67 | ("R1 (A)","R2 (A)","R3 (A)","R4 (A)","R5 (A)","R6 (A)","rho1 (cm-2)","rho2 (cm-2)","rho3 (cm-2)","rho4 (cm-2)","rho5 (cm-2)","rho6 (cm-2)","rho7 (cm-2)","L (A)","Rhos (cm-3)"), |
---|
68 | ("%1.3e","%1.3e","%1.3e","%1.3e","%1.3e","%1.3e","%1.3e","%1.3e","%1.3e","%1.3e","%1.3e","%1.3e","%1.3e","%1.3e","%1.3e"), |
---|
69 | (True,True,True,True,True,True,True,True,True,True,True,True,True,True,True))] |
---|
70 | |
---|
71 | from LSsca |
---|
72 | ''' |
---|
73 | def __init__(self): |
---|
74 | Model.__init__(self) |
---|
75 | self.IntensityFunc=self.PcylmultiminFunction #function |
---|
76 | self.N=0 |
---|
77 | self.q=Qlogspace(1e-4,1.,500.) #q range(x scale) |
---|
78 | self.Arg=[13.,14.,15.,16.,17.,18.,9.418e+010,8.798e+010,4.041e+011,3.925e+011,4.224e+011,3.902e+011,9.418e+010,400.,1e10] #list of defaults parameters |
---|
79 | self.Format=["%1.3e","%1.3e","%1.3e","%1.3e","%1.3e","%1.3e","%1.3e","%1.3e","%1.3e","%1.3e","%1.3e","%1.3e","%1.3e","%1.3e","%1.3e"] #list of c format |
---|
80 | self.istofit=[True,True,True,True,True,True,True,True,True,True,True,True,True,True,True] #list of boolean for fitting |
---|
81 | self.name="Cylinder with six levels" #name of the model |
---|
82 | self.Doc=["R1 (A)","R2 (A)","R3 (A)","R4 (A)","R5 (A)","R6 (A)","rho1 (cm-2)","rho2 (cm-2)","rho3 (cm-2)","rho4 (cm-2)","rho5 (cm-2)","rho6 (cm-2)","rho7 (cm-2)","L (A)","Rhos (cm-3)"] #list of description for parameters |
---|
83 | |
---|
84 | if __name__=="__main__": |
---|
85 | ''' |
---|
86 | test code |
---|
87 | ''' |
---|
88 | modl=Pcylmultimin() |
---|
89 | #plot the model |
---|
90 | import Gnuplot |
---|
91 | gp=Gnuplot.Gnuplot() |
---|
92 | gp("set logscale xy") |
---|
93 | c=Gnuplot.Data(modl.q,modl.getIntensity(),with_='points') |
---|
94 | gp.plot(c) |
---|
95 | raw_input("enter") |
---|
96 | #plot and fit the noisy model |
---|
97 | yn=modl.getNoisy(0.4) |
---|
98 | cn=Gnuplot.Data(modl.q,yn,with_='points') |
---|
99 | res=modl.fit(yn) |
---|
100 | cf=Gnuplot.Data(modl.q,modl.IntensityFunc(modl.q,res),with_='lines') |
---|
101 | gp.plot(c,cn,cf) |
---|
102 | raw_input("enter") |
---|
103 | #plot and fit the noisy model with fitBounds |
---|
104 | bounds=modl.getBoundsFromParam() #[250.0,2e11,1e10,1.5e15] |
---|
105 | res2=modl.fitBounds(yn,bounds) |
---|
106 | print res2 |
---|
107 | raw_input("enter") |
---|