Ticket #678: RandomToOrientedCoreShellShellChains.py

File RandomToOrientedCoreShellShellChains.py, 15.1 KB (added by pkienzle, 8 years ago)
Line 
1"""
2Example model
3"""
4
5from math import *
6import os
7import sys
8import numpy
9
10#name
11name = "RandomToOrientedCoreShellShellChains"
12
13#title
14title = "RandomToOrientedCoreShellShellChains"
15
16#description
17description = "Core-Shell-Shell Linear Chains Random-To-Oriented About X"
18
19#parameters
20parameters = [ 
21                ['core_sld', 'Inv. Ang^2', 6.9E-6, [-numpy.inf, numpy.inf], '', ''],
22                ['magcore_sld', 'Inv. Ang^2', 1.5E-6, [-numpy.inf, numpy.inf], '', ''],
23                ['shell1_sld', 'Inv. Ang^2', 1.0E-7, [-numpy.inf, numpy.inf], '', ''],
24                ['shell2_sld', 'Inv. Ang^2', 1.0E-7, [-numpy.inf, numpy.inf], '', ''],
25                ['solvent_sld', 'Inv. Ang^2', 6.35E-6, [-numpy.inf, numpy.inf], '', ''],
26                ['volume_fraction', '', 0.0004, [0.0, numpy.inf], '', ''],
27                ['normalization_radius', '', 125, [0.0, numpy.inf], '', ''],
28                ['singlet_fraction', '', 0.5, [0.0, numpy.inf], '', ''],
29                ['dimer_fraction', '', 0.3, [0.0, numpy.inf], '', ''],
30                ['trimer_fraction', '', 0.2, [0.0, numpy.inf], '', ''],
31                ['quadramer_fraction', '', 0.1, [0.0, numpy.inf], '', ''],
32                ['pentamer_fraction', '', 0.1, [0.0, numpy.inf], '', ''],
33                ['core_radius', 'Ang.', 125, [0.0, numpy.inf], 'volume', ''],
34                ['shell1_thickness', 'Ang.', 25, [0.0, numpy.inf], 'volume', ''],
35                ['shell2_thickness', 'Ang.', 25, [0.0, numpy.inf], 'volume', ''],
36                ['length', 'Sphere to sphere length in Ang.', 350, [0.0, numpy.inf], 'volume', ''],
37                ['m_var', '1=alongchain;2=alongchainskewtoX,3=alongfield', 1, [0, 3], '', ''],
38                ['viewing_angle', 'w.r.t. applied field (degrees)', 0, [-90, 90], '', ''],
39                ['sigma', 'Standard deviation of chain orientation (degrees)', 60, [0.1, 5000], '', ''],
40                ['cross_section', '0=Unpol,1=UU,2=UD,3=DD,4=DU,5=DD-UU,6=|DD-UU|,7=DU+UD,8=DU-UD,9=|DU-UD|', 0.0, [0, 10], '', ''],
41                ['porod_slope', '', -4.0, [-numpy.inf, numpy.inf], '', ''],
42                ['porod_scale', '', 0.0, [0.0, numpy.inf], '', '']
43             ]
44
45form_volume = """
46    return 1.0;
47"""
48
49def ER(*arg): 
50    return 1.0 
51
52Iq = """
53
54    double Vol = M_4PI_3*pow(normalization_radius,3);
55    if(Vol == 0){Vol = 1E-10;}
56    double Q_X = q*cos(viewing_angle*M_PI_180);
57    double Q_Y = q*sin(viewing_angle*M_PI_180);
58
59    double GammaMatrix[3][3];
60    double QA_hat= cos(viewing_angle*M_PI_180);
61    double QB_hat= sin(viewing_angle*M_PI_180);
62    double QC_hat= 0.0;
63    GammaMatrix[0][0] = 1.0 - QA_hat*QA_hat; GammaMatrix[0][1] = -QA_hat*QB_hat; GammaMatrix[0][2] = -QA_hat*QC_hat;
64    GammaMatrix[1][0] = 1.0 - QB_hat*QB_hat; GammaMatrix[1][1] = -QA_hat*QB_hat; GammaMatrix[1][2] = -QB_hat*QC_hat;
65    GammaMatrix[2][0] = 1.0 - QC_hat*QC_hat; GammaMatrix[2][1] = -QA_hat*QC_hat; GammaMatrix[2][2] = -QB_hat*QC_hat;
66
67
68    double AmpR1 = (sin(q*core_radius)/q-core_radius*cos(q*core_radius))/(q*q);
69    double AmpR2 = (sin(q*(core_radius + shell1_thickness))/q-(core_radius + shell1_thickness)*cos(q*(core_radius + shell1_thickness)))/(q*q);
70    double AmpR3 = (sin(q*(core_radius + shell1_thickness + shell2_thickness))/q-(core_radius + shell1_thickness + shell2_thickness)*cos(q*(core_radius + shell1_thickness + shell2_thickness)))/(q*q);
71    double Amp = 4.0*M_PI*((core_sld-solvent_sld)*AmpR1 + (shell1_sld-solvent_sld)*(AmpR2-AmpR1) + (shell2_sld-solvent_sld)*(AmpR3-AmpR2));
72    double MAmp = 4.0*M_PI*(magcore_sld)*AmpR1;
73
74    const int a_steps = 17;
75    const int b_steps = 12;
76    double ChainProjX[20][20];//must be larger than a_steps, b_steps
77    double ChainProjY[20][20];//must be larger than a_steps, b_steps
78    double ChainProjZ[20][20];//must be larger than a_steps, b_steps
79    double Norm = 0.0;
80    double quad_place = 0.0;
81    double angle_from_VA, angle_from_x, variable1, variable2;//angle is defined relative to viewing_angle (not x-axis)
82    for(int a=0; a<a_steps; a++){//(a*2.0+1.0-90.0)*M_PI_180 for 0->90
83        for(int b=0; b<b_steps; b++){//(b*45.0)*M_PI_180 for 0-> 12
84            angle_from_VA = (a*5.0-90.0)*M_PI_180;
85            double phi = (b*45.0)*M_PI_180;
86
87                //********************************************
88                //Rotate chains:
89                double alpha[3], beta[3], gamma[3];
90                alpha[0] = cos(viewing_angle*M_PI_180); alpha[1] = sin(viewing_angle*M_PI_180); alpha[2] = 0;
91                beta[0] = cos(viewing_angle*M_PI_180+M_PI/2); beta[1] = sin(viewing_angle*M_PI_180+M_PI/2); beta[2] = 0;
92                gamma[0] = 0.0; gamma[1] = 0.0; gamma[2] = 1.0;
93                int points = 3;
94                double LHS_Matrix[3][3];
95                LHS_Matrix[0][0] = alpha[0]; LHS_Matrix[0][1] = alpha[1]; LHS_Matrix[0][2] = alpha[2];
96                LHS_Matrix[1][0] = beta[0]; LHS_Matrix[1][1] = beta[1]; LHS_Matrix[1][2] = beta[2];
97                LHS_Matrix[2][0] = gamma[0]; LHS_Matrix[2][1] = gamma[1]; LHS_Matrix[2][2] = gamma[2];
98                double RHS_Matrix[3] = {cos(angle_from_VA), sin(angle_from_VA)*cos(phi), sin(angle_from_VA)*sin(phi)};
99                double Alignment[3];
100
101                if(alpha[0] <= 0.01 && alpha[0] >= -0.01){
102                LHS_Matrix[1][0] = alpha[0]; LHS_Matrix[1][1] = alpha[1]; LHS_Matrix[1][2] = alpha[2];
103                LHS_Matrix[0][0] = beta[0]; LHS_Matrix[0][1] = beta[1]; LHS_Matrix[0][2] = beta[2];
104                RHS_Matrix[0] = sin(angle_from_VA)*cos(phi);
105                RHS_Matrix[1] = cos(angle_from_VA);
106                }
107               
108
109                //First row Gauss-Jordan ellimination:
110                for(int col = 0; col<points-1; col++){
111                for(int row = col+1; row<points; row++){
112                for(int p=1; p<points-col; p++){
113                        LHS_Matrix[row][col+p] = LHS_Matrix[row][col+p] -(LHS_Matrix[row][col] / LHS_Matrix[col][col])*LHS_Matrix[col][col+p];
114                }//end p loop
115                        RHS_Matrix[row] = RHS_Matrix[row] -(LHS_Matrix[row][col] / LHS_Matrix[col][col])*RHS_Matrix[col];
116                        LHS_Matrix[row][col] = LHS_Matrix[row][col] -(LHS_Matrix[row][col] / LHS_Matrix[col][col])*LHS_Matrix[col][col];
117                }//end row loop
118                }//end col loop
119
120                //Additional rows of Gauss-Jordan ellimination:
121                for(int col = 1; col<points; col++){
122                for(int row = 0; row<col; row++){
123                for(int p=1; p<points-col; p++){
124                        LHS_Matrix[row][col+p] = LHS_Matrix[row][col+p] -(LHS_Matrix[row][col] / LHS_Matrix[col][col])*LHS_Matrix[col][col+p];
125                }//end p loop
126                        RHS_Matrix[row] = RHS_Matrix[row] -(LHS_Matrix[row][col] / LHS_Matrix[col][col])*RHS_Matrix[col];
127                        LHS_Matrix[row][col] = LHS_Matrix[row][col] -(LHS_Matrix[row][col] / LHS_Matrix[col][col])*LHS_Matrix[col][col];
128                }//end row loop
129                }//end col loop
130
131
132                for(int p=0; p<points; p++){
133                Alignment[p] = RHS_Matrix[p] / LHS_Matrix[p][p];
134                }//end p loop
135
136                double Unit = sqrt(Alignment[0]*Alignment[0] + Alignment[1]*Alignment[1] + Alignment[2]*Alignment[2]);
137                ChainProjX[a][b] = Alignment[0]/Unit;
138                ChainProjY[a][b] = Alignment[1]/Unit;
139                ChainProjZ[a][b] = Alignment[2]/Unit;
140                //********************************************
141
142            //angle_from_x = viewing_angle + angle_from_VA;
143            angle_from_x = acos(ChainProjX[a][b]);
144            if(angle_from_x >= -M_PI && angle_from_x <= M_PI){quad_place = 0.0;}
145            if(angle_from_x < -M_PI){quad_place = M_PI;}
146            if(angle_from_x > M_PI){quad_place = -M_PI;}
147            variable1 = (angle_from_x + quad_place - 0.0)/(sigma*M_PI_180);
148            variable2 = sqrt(2.0*M_PI)*sigma*180.0/M_PI;
149            Norm += sqrt(pow(sin(angle_from_VA),2))*(exp( -0.5*pow( variable1,2) ) / variable2);
150        }
151    }
152
153    double CrossSection[4][5] = {0.0};
154    double Nucl = 0.0;
155    double Nucl2 = 0.0;
156
157    double anglewt;
158    for(int a=0; a<17; a++){
159        for(int b=0; b<12; b++){
160            angle_from_VA = (a*5.0-90.0)*M_PI_180;
161            double phi = (b*45.0)*M_PI_180;
162            //angle_from_x = viewing_angle + angle_from_VA;
163            angle_from_x = acos(ChainProjX[a][b]);
164            if(angle_from_x >= -M_PI && angle_from_x <= M_PI){quad_place = 0.0;}
165            if(angle_from_x < -M_PI){quad_place = M_PI;}
166            if(angle_from_x > M_PI){quad_place = -M_PI;}
167            variable1 = (angle_from_x + quad_place - 0.0)/(sigma*M_PI_180);
168            variable2 = sqrt(2.0*M_PI)*sigma/M_PI_180;
169            anglewt = sqrt(pow(sin(angle_from_VA),2))*(exp( -0.5*pow( variable1,2) ) / variable2)/Norm;
170
171            double mag_a_real = 0.0;
172            double mag_a_img = 0.0;
173            double mag_b_real = 0.0;
174            double mag_b_img = 0.0;
175            double mag_c_real = 0.0;
176            double mag_c_img = 0.0;
177            double nucl_real = 0.0;
178            double nucl_img = 0.0;
179            double A_real = 0.0;
180            double A_img = 0.0;
181            double B_real = 0.0;
182            double B_img = 0.0;
183            double C_real = 0.0;
184            double C_img = 0.0;
185
186            for(int k=0; k<5; k++){
187                double number = 1.0*k + 1.0;
188                nucl_real += Amp*cos(k*length*(Q_X*ChainProjX[a][b] + Q_Y*ChainProjY[a][b]));
189                nucl_img += Amp*sin(k*length*(Q_X*ChainProjX[a][b] + Q_Y*ChainProjY[a][b]));
190                if(m_var <= 1.5){
191                mag_a_real = MAmp*ChainProjX[a][b]*cos(k*length*(Q_X*ChainProjX[a][b] + Q_Y*ChainProjY[a][b]));
192                mag_b_real = MAmp*ChainProjY[a][b]*cos(k*length*(Q_X*ChainProjX[a][b] + Q_Y*ChainProjY[a][b]));
193                mag_c_real = MAmp*ChainProjZ[a][b]*cos(k*length*(Q_X*ChainProjX[a][b] + Q_Y*ChainProjY[a][b]));
194                mag_a_img = MAmp*ChainProjX[a][b]*sin(k*length*(Q_X*ChainProjX[a][b] + Q_Y*ChainProjY[a][b]));
195                mag_b_img = MAmp*ChainProjY[a][b]*sin(k*length*(Q_X*ChainProjX[a][b] + Q_Y*ChainProjY[a][b]));
196                mag_c_img = MAmp*ChainProjZ[a][b]*sin(k*length*(Q_X*ChainProjX[a][b] + Q_Y*ChainProjY[a][b]));
197                }
198                if(m_var > 1.5 && m_var <= 2.5){
199                mag_a_real = MAmp*sqrt(ChainProjX[a][b]*ChainProjX[a][b])*cos(k*length*(Q_X*ChainProjX[a][b] + Q_Y*ChainProjY[a][b]));
200                mag_b_real = MAmp*ChainProjY[a][b]*cos(k*length*(Q_X*ChainProjX[a][b] + Q_Y*ChainProjY[a][b]));
201                mag_c_real = MAmp*ChainProjZ[a][b]*cos(k*length*(Q_X*ChainProjX[a][b] + Q_Y*ChainProjY[a][b]));
202                mag_a_img = MAmp*sqrt(ChainProjX[a][b]*ChainProjX[a][b])*sin(k*length*(Q_X*ChainProjX[a][b] + Q_Y*ChainProjY[a][b]));
203                mag_b_img = MAmp*ChainProjY[a][b]*sin(k*length*(Q_X*ChainProjX[a][b] + Q_Y*ChainProjY[a][b]));
204                mag_c_img = MAmp*ChainProjZ[a][b]*sin(k*length*(Q_X*ChainProjX[a][b] + Q_Y*ChainProjY[a][b]));
205                }
206                if(m_var > 2.5 && m_var <= 3.5){
207                mag_a_real = MAmp*cos(k*length*(Q_X*ChainProjX[a][b] + Q_Y*ChainProjY[a][b]));
208                mag_b_real = 0.0;
209                mag_c_real = 0.0;
210                mag_a_img = MAmp*sin(k*length*(Q_X*ChainProjX[a][b] + Q_Y*ChainProjY[a][b]));
211                mag_b_img = 0.0;
212                mag_c_img = 0.0;
213                }
214                A_real = GammaMatrix[0][0]*mag_a_real + GammaMatrix[0][1]*mag_b_real + GammaMatrix[0][2]*mag_c_real;
215                A_img = GammaMatrix[0][0]*mag_a_img + GammaMatrix[0][1]*mag_b_img + GammaMatrix[0][2]*mag_c_img;
216                B_real = GammaMatrix[1][0]*mag_b_real + GammaMatrix[1][1]*mag_a_real + GammaMatrix[1][2]*mag_c_real;
217                B_img = GammaMatrix[1][0]*mag_b_img + GammaMatrix[1][1]*mag_a_img + GammaMatrix[1][2]*mag_c_img;
218                C_real = GammaMatrix[2][0]*mag_c_real + GammaMatrix[2][1]*mag_a_real + GammaMatrix[2][2]*mag_b_real;
219                C_img = GammaMatrix[2][0]*mag_c_img + GammaMatrix[2][1]*mag_a_img + GammaMatrix[2][2]*mag_b_img;
220                CrossSection[0][k] += 0.5*anglewt*(pow(nucl_real-A_real,2)+pow(nucl_img-A_img,2))/(Vol*number);
221                CrossSection[1][k] += 0.5*anglewt*(pow(-B_real+C_img,2)+pow(-B_img-C_real,2))/(Vol*number);
222                CrossSection[2][k] += 0.5*anglewt*(pow(nucl_real+A_real,2)+pow(nucl_img+A_img,2))/(Vol*number);
223                CrossSection[3][k] += 0.5*anglewt*(pow(-B_real-C_img,2)+pow(-B_img+C_real,2))/(Vol*number);
224                Nucl2 += anglewt*(pow(nucl_real-A_real,2)+pow(nucl_img-A_img,2))/(Vol*number);
225                Nucl += anglewt*(pow(nucl_real,2)+pow(nucl_img,2))/(Vol*number);
226
227            }
228       }
229    }
230
231    double Unpol = 0.0;
232    double Intensity = 0.0;
233    double FractionScale = singlet_fraction + dimer_fraction + trimer_fraction + quadramer_fraction + pentamer_fraction;
234    if(FractionScale == 0){FractionScale = 1.0;}
235
236    for(int n=0; n<4; n++){
237    Unpol += singlet_fraction*CrossSection[n][0] + dimer_fraction*CrossSection[n][1] + trimer_fraction*CrossSection[n][2] + quadramer_fraction*CrossSection[n][3] + pentamer_fraction*CrossSection[n][4];
238    }
239
240    if(cross_section < 0.5){
241    Intensity = Unpol;
242    }
243
244    if(cross_section >= 0.5 && cross_section <= 1.5){
245    int n = 0;
246    Intensity = singlet_fraction*CrossSection[n][0] + dimer_fraction*CrossSection[n][1] + trimer_fraction*CrossSection[n][2] + quadramer_fraction*CrossSection[n][3] + pentamer_fraction*CrossSection[n][4];
247    }
248
249    if(cross_section >= 1.5 && cross_section <= 2.5){
250    int n = 1;
251    Intensity = singlet_fraction*CrossSection[n][0] + dimer_fraction*CrossSection[n][1] + trimer_fraction*CrossSection[n][2] + quadramer_fraction*CrossSection[n][3] + pentamer_fraction*CrossSection[n][4];
252    }
253
254    if(cross_section >= 2.5 && cross_section <= 3.5){
255    int n = 2;
256    Intensity = singlet_fraction*CrossSection[n][0] + dimer_fraction*CrossSection[n][1] + trimer_fraction*CrossSection[n][2] + quadramer_fraction*CrossSection[n][3] + pentamer_fraction*CrossSection[n][4];
257    }
258
259    if(cross_section >= 3.5 && cross_section <= 4.5){
260    int n = 3;
261    Intensity = singlet_fraction*CrossSection[n][0] + dimer_fraction*CrossSection[n][1] + trimer_fraction*CrossSection[n][2] + quadramer_fraction*CrossSection[n][3] + pentamer_fraction*CrossSection[n][4];
262    }
263
264    if(cross_section >= 4.5 && cross_section <= 6.5){
265    int n = 2;
266    Intensity = singlet_fraction*CrossSection[n][0] + dimer_fraction*CrossSection[n][1] + trimer_fraction*CrossSection[n][2] + quadramer_fraction*CrossSection[n][3] + pentamer_fraction*CrossSection[n][4];
267    n = 0;
268    Intensity = Intensity - (singlet_fraction*CrossSection[n][0] + dimer_fraction*CrossSection[n][1] + trimer_fraction*CrossSection[n][2] + quadramer_fraction*CrossSection[n][3] + pentamer_fraction*CrossSection[n][4]);
269        if(cross_section >= 5.5){double holder = pow(Intensity,2); Intensity = sqrt(holder);}
270    }
271
272    if(cross_section >= 6.5 && cross_section < 7.5){
273    int n = 3;
274    Intensity = singlet_fraction*CrossSection[n][0] + dimer_fraction*CrossSection[n][1] + trimer_fraction*CrossSection[n][2] + quadramer_fraction*CrossSection[n][3] + pentamer_fraction*CrossSection[n][4];
275    n = 1;
276    Intensity = Intensity + (singlet_fraction*CrossSection[n][0] + dimer_fraction*CrossSection[n][1] + trimer_fraction*CrossSection[n][2] + quadramer_fraction*CrossSection[n][3] + pentamer_fraction*CrossSection[n][4]);
277    }
278
279    if(cross_section >= 7.5){
280    int n = 3;
281    Intensity = singlet_fraction*CrossSection[n][0] + dimer_fraction*CrossSection[n][1] + trimer_fraction*CrossSection[n][2] + quadramer_fraction*CrossSection[n][3] + pentamer_fraction*CrossSection[n][4];
282    n = 1;
283    Intensity = Intensity - (singlet_fraction*CrossSection[n][0] + dimer_fraction*CrossSection[n][1] + trimer_fraction*CrossSection[n][2] + quadramer_fraction*CrossSection[n][3] + pentamer_fraction*CrossSection[n][4]);
284        if(cross_section >= 8.5){double holder = pow(Intensity,2); Intensity = sqrt(holder);}
285    }
286
287    if(Intensity*1E6 < Unpol){
288    Intensity = 0.0;
289    }
290
291    double PorodSlope = porod_scale*pow(q,porod_slope);
292
293    return Intensity*volume_fraction*(1E8)/FractionScale + PorodSlope;
294
295"""
296
297Iqxy = """
298    return Iq(sqrt(qx*qx+qy*qy), core_sld, magcore_sld, shell1_sld, shell2_sld, solvent_sld, volume_fraction, normalization_radius, singlet_fraction, dimer_fraction, trimer_fraction, quadramer_fraction, pentamer_fraction, core_radius, shell1_thickness, shell2_thickness, length, m_var, viewing_angle, sigma, cross_section, porod_slope, porod_scale);
299"""