source: sasmodels/sasmodels/models/lamellar_hg_stack_caille.c @ d86f0fc

core_shell_microgelsmagnetic_modelticket-1257-vesicle-productticket_1156ticket_1265_superballticket_822_more_unit_tests
Last change on this file since d86f0fc was 1c6286d, checked in by Paul Kienzle <pkienzle@…>, 7 years ago

lamellar stack/hg-stack caille: code tidying, changing loop terminator from 'i⇐N-1' to the more usual 'i<N'

  • Property mode set to 100644
File size: 1.7 KB
Line 
1/*      LamellarCailleHG kernel - allows for name changes of passed parameters ...
2    Maths identical to LamellarCaille apart from the line for P(Q)
3*/
4
5static double
6Iq(double qval,
7   double length_tail,
8   double length_head,
9   double fp_Nlayers,
10   double dd,
11   double Cp,
12   double tail_sld,
13   double head_sld,
14   double solvent_sld)
15{
16  int Nlayers = (int)(fp_Nlayers+0.5);    //cast to an integer for the loop
17  double inten,Pq,Sq,alpha,temp,t2;
18  //double dQ, dQDefault, t1, t3;
19  // from wikipedia 0.577215664901532860606512090082402431042159335
20  const double Euler = 0.577215664901533;   // Euler's constant, increased sig figs for new models Feb 2015
21  //dQDefault = 0.0;    //[=] 1/A, q-resolution, default value
22  //dQ = dQDefault; // REMOVED UNUSED dQ calculations for new models Feb 2015
23
24  Pq = (head_sld-solvent_sld)*(sin(qval*(length_head+length_tail))-sin(qval*length_tail))
25       + (tail_sld-solvent_sld)*sin(qval*length_tail);
26  Pq *= Pq;
27  Pq *= 4.0/(qval*qval);
28
29  Sq = 0.0;
30  for(int ii=1; ii < Nlayers; ii++) {
31    temp = 0.0;
32    alpha = Cp/4.0/M_PI/M_PI*(log(M_PI*ii) + Euler);
33    //t1 = 2.0*dQ*dQ*dd*dd*alpha;
34    t2 = 2.0*qval*qval*dd*dd*alpha;
35    //t3 = dQ*dQ*dd*dd*ii*ii;
36
37    temp = 1.0-(double)ii/(double)Nlayers;
38    //temp *= cos(dd*qval*ii/(1.0+t1));
39    temp *= cos(dd*qval*ii);
40    //if (temp < 0) printf("q=%g: ii=%d, cos(dd*q*ii)=cos(%g) < 0\n",qval,ii,dd*qval*ii);
41    //temp *= exp(-1.0*(t2 + t3)/(2.0*(1.0+t1)) );
42    temp *= exp(-t2/2.0);
43    //temp /= sqrt(1.0+t1);
44
45    Sq += temp;
46  }
47
48  Sq *= 2.0;
49  Sq += 1.0;
50
51  //if (Sq < 0) printf("q=%g: S(q) =%g\n", qval, Sq);
52
53  inten = 2.0*M_PI*Pq*Sq/(dd*qval*qval);
54
55  inten *= 1.0e-04;   // 1/A to 1/cm
56  return inten;
57}
58
Note: See TracBrowser for help on using the repository browser.