# Ticket #1080: PCylinder_six_layers - Copy.py

File PCylinder_six_layers - Copy.py, 5.2 KB (added by smk78, 6 years ago)

pySAXS multi-shell cylinder model

Line
1from model import Model
2from pySAXS.LS.LSsca import Qlogspace
3from pySAXS.LS.LSsca import *
4import numpy
5
6class 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
84if __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")