Changeset 4329539 in sasmodels


Ignore:
Timestamp:
Mar 5, 2017 4:23:38 PM (7 years ago)
Author:
butler
Branches:
master, core_shell_microgels, costrafo411, magnetic_model, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
bb73096
Parents:
e1aa129
Message:

add comments to rpa c code to help code readability. Addresses ticket
#593

File:
1 edited

Legend:

Unmodified
Added
Removed
  • sasmodels/models/rpa.c

    r6351bfa r4329539  
    3939  double S31,S32,S33,S34,S41,S42,S43,S44; 
    4040  double Lad,Lbd,Lcd,Nav,Intg; 
    41  
     41   
     42  // Set values for non existent parameters (eg. no A or B in case 0 and 1 etc) 
    4243  //icase was shifted to N-1 from the original code 
    4344  if (icase <= 1){ 
     
    5758    Kab = Kac = Kad = -0.0004; 
    5859  } 
    59  
     60  
     61  //set up values for cross terms in case of block copolymers (1,3,4,6,7,8,9) 
    6062  Nab=sqrt(N[0]*N[1]); 
    6163  Nac=sqrt(N[0]*N[2]); 
     
    7981  Phicd=sqrt(Phi[2]*Phi[3]); 
    8082 
     83  // Calculate Q^2 * Rg^2 for each homopolymer assuming random walk  
    8184  Xa=q*q*b[0]*b[0]*N[0]/6.0; 
    8285  Xb=q*q*b[1]*b[1]*N[1]/6.0; 
     
    8487  Xd=q*q*b[3]*b[3]*N[3]/6.0; 
    8588 
    86   Paa=2.0*(exp(-Xa)-1.0+Xa)/(Xa*Xa); 
    87   S0aa=N[0]*Phi[0]*v[0]*Paa; 
    88   Pab=((1.0-exp(-Xa))/Xa)*((1.0-exp(-Xb))/Xb); 
     89  //calculate all partial structure factors Pij and normalize n^2 
     90  Paa=2.0*(exp(-Xa)-1.0+Xa)/(Xa*Xa); // free A chain form factor 
     91  S0aa=N[0]*Phi[0]*v[0]*Paa; // Phi * Vp * P(Q)= I(Q0)/delRho^2 
     92  Pab=((1.0-exp(-Xa))/Xa)*((1.0-exp(-Xb))/Xb); //AB diblock (anchored Paa * anchored Pbb) partial form factor 
    8993  S0ab=(Phiab*vab*Nab)*Pab; 
    90   Pac=((1.0-exp(-Xa))/Xa)*exp(-Xb)*((1.0-exp(-Xc))/Xc); 
     94  Pac=((1.0-exp(-Xa))/Xa)*exp(-Xb)*((1.0-exp(-Xc))/Xc); //ABC triblock AC partial form factor 
    9195  S0ac=(Phiac*vac*Nac)*Pac; 
    92   Pad=((1.0-exp(-Xa))/Xa)*exp(-Xb-Xc)*((1.0-exp(-Xd))/Xd); 
     96  Pad=((1.0-exp(-Xa))/Xa)*exp(-Xb-Xc)*((1.0-exp(-Xd))/Xd); //ABCD four block 
    9397  S0ad=(Phiad*vad*Nad)*Pad; 
    9498 
    9599  S0ba=S0ab; 
    96   Pbb=2.0*(exp(-Xb)-1.0+Xb)/(Xb*Xb); 
     100  Pbb=2.0*(exp(-Xb)-1.0+Xb)/(Xb*Xb); // free B chain 
    97101  S0bb=N[1]*Phi[1]*v[1]*Pbb; 
    98   Pbc=((1.0-exp(-Xb))/Xb)*((1.0-exp(-Xc))/Xc); 
     102  Pbc=((1.0-exp(-Xb))/Xb)*((1.0-exp(-Xc))/Xc); // BC diblock 
    99103  S0bc=(Phibc*vbc*Nbc)*Pbc; 
    100   Pbd=((1.0-exp(-Xb))/Xb)*exp(-Xc)*((1.0-exp(-Xd))/Xd); 
     104  Pbd=((1.0-exp(-Xb))/Xb)*exp(-Xc)*((1.0-exp(-Xd))/Xd); // BCD triblock 
    101105  S0bd=(Phibd*vbd*Nbd)*Pbd; 
    102106 
    103107  S0ca=S0ac; 
    104108  S0cb=S0bc; 
    105   Pcc=2.0*(exp(-Xc)-1.0+Xc)/(Xc*Xc); 
     109  Pcc=2.0*(exp(-Xc)-1.0+Xc)/(Xc*Xc); // Free C chain 
    106110  S0cc=N[2]*Phi[2]*v[2]*Pcc; 
    107   Pcd=((1.0-exp(-Xc))/Xc)*((1.0-exp(-Xd))/Xd); 
     111  Pcd=((1.0-exp(-Xc))/Xc)*((1.0-exp(-Xd))/Xd); // CD diblock 
    108112  S0cd=(Phicd*vcd*Ncd)*Pcd; 
    109113 
     
    111115  S0db=S0bd; 
    112116  S0dc=S0cd; 
    113   Pdd=2.0*(exp(-Xd)-1.0+Xd)/(Xd*Xd); 
     117  Pdd=2.0*(exp(-Xd)-1.0+Xd)/(Xd*Xd); // free D chain 
    114118  S0dd=N[3]*Phi[3]*v[3]*Pdd; 
     119 
     120  // Reset all unused partial structure factors to 0 (depends on case) 
    115121  //icase was shifted to N-1 from the original code 
    116122  switch(icase){ 
     
    193199  S0dc=S0cd; 
    194200 
     201  // self chi parameter is 0 ... of course 
    195202  Kaa=0.0; 
    196203  Kbb=0.0; 
     
    243250  ZZ=S0ad*(T11*S0ad+T12*S0bd+T13*S0cd)+S0bd*(T21*S0ad+T22*S0bd+T23*S0cd)+S0cd*(T31*S0ad+T32*S0bd+T33*S0cd); 
    244251 
     252  // D is considered the matrix or background component so enters here 
    245253  m=1.0/(S0dd-ZZ); 
    246254 
     
    297305  S44=S11+S22+S33+2.0*S12+2.0*S13+2.0*S23; 
    298306 
     307  //calculate contrast where L[i] is the scattering length of i and D is the matrix 
     308  //need to verify why the sqrt of Nav rather than just Nav (assuming v is molar volume) 
    299309  Nav=6.022045e+23; 
    300310  Lad=(L[0]/v[0]-L[3]/v[3])*sqrt(Nav); 
Note: See TracChangeset for help on using the changeset viewer.