Changeset 32e3c9b in sasmodels for sasmodels/kernel_iq.c


Ignore:
Timestamp:
Jul 21, 2016 2:08:04 PM (8 years ago)
Author:
Paul Kienzle <pkienzle@…>
Branches:
master, core_shell_microgels, costrafo411, magnetic_model, release_v0.94, release_v0.95, ticket-1257-vesicle-product, ticket_1156, ticket_1265_superball, ticket_822_more_unit_tests
Children:
b966a96
Parents:
42356c8
Message:

dll version of magnetic sld

File:
1 edited

Legend:

Unmodified
Added
Removed

  • sasmodels/kernel_iq.c

    ra738209 r32e3c9b  
    3333#endif 
    3434 
     35#ifdef MAGNETIC 
     36const int32_t magnetic[] = { MAGNETIC_PARS }; 
     37#endif 
     38 
     39#ifdef MAGNETIC 
     40// Return value restricted between low and high 
     41static double clip(double value, double low, double high) 
     42{ 
     43    return (value < low ? low : (value > high ? high : value)); 
     44} 
     45 
     46// Compute spin cross sections given in_spin and out_spin 
     47// To convert spin cross sections to sld b: 
     48//     uu * (sld - m_sigma_x); 
     49//     dd * (sld + m_sigma_x); 
     50//     ud * (m_sigma_y + 1j*m_sigma_z); 
     51//     du * (m_sigma_y - 1j*m_sigma_z); 
     52static void spins(double in_spin, double out_spin, 
     53    double *uu, double *dd, double *ud, double *du) 
     54{ 
     55    in_spin = clip(in_spin, 0.0, 1.0); 
     56    out_spin = clip(out_spin, 0.0, 1.0); 
     57        *uu = sqrt(sqrt(in_spin * out_spin)); 
     58        *dd = sqrt(sqrt((1.0-in_spin) * (1.0-out_spin))); 
     59        *ud = sqrt(sqrt(in_spin * (1.0-out_spin))); 
     60        *du = sqrt(sqrt((1.0-in_spin) * out_spin)); 
     61} 
     62 
     63// Convert polar to rectangular coordinates. 
     64static void polrec(double r, double theta, double phi, 
     65    double *x, double *y, double *z) 
     66{ 
     67    double cos_theta, sin_theta, cos_phi, sin_phi; 
     68    SINCOS(theta*M_PI_180, sin_theta, cos_theta); 
     69    SINCOS(phi*M_PI_180, sin_phi, cos_phi); 
     70    *x = r * cos_theta * cos_phi; 
     71    *y = r * sin_theta; 
     72    *z = -r * cos_theta * sin_phi; 
     73} 
     74#endif 
    3575 
    3676kernel 
    … …  
    4080    const int32_t pd_stop,      // where we are stopping in the polydispersity loop 
    4181    global const ProblemDetails *details, 
    42     global const double *values, 
     82   // global const  // TODO: make it const again! 
     83    double *values, 
    4384    global const double *q, // nq q values, with padding to boundary 
    4485    global double *result,  // nq+3 return values, again with padding 
    … …  
    5293 
    5394  // Fill in the initial variables 
     95  //   values[0] is scale 
     96  //   values[1] is background 
    5497  #ifdef USE_OPENMP 
    5598  #pragma omp parallel for 
    … …  
    58101    pvec[k] = values[k+2]; 
    59102  } 
     103#ifdef MAGNETIC 
     104  const double up_frac_i = values[NPARS+2]; 
     105  const double up_frac_f = values[NPARS+3]; 
     106  const double up_angle = values[NPARS+4]; 
     107  #define MX(_k) (values[NPARS+5+3*_k]) 
     108  #define MY(_k) (values[NPARS+6+3*_k]) 
     109  #define MZ(_k) (values[NPARS+7+3*_k]) 
     110 
     111  // TODO: precompute this on the python side 
     112  // Convert polar to rectangular coordinates in place. 
     113  if (pd_start == 0) {  // Update in place; only do this for the first hunk! 
     114//printf("spin: %g %g %g\n", up_frac_i, up_frac_f, up_angle); 
     115    for (int mag=0; mag < NUM_MAGNETIC; mag++) { 
     116//printf("mag %d: %g %g %g\n", mag, MX(mag), MY(mag), MZ(mag)); 
     117        polrec(MX(mag), MY(mag), MZ(mag), &MX(mag), &MY(mag), &MZ(mag)); 
     118//printf("   ==>: %g %g %g\n", MX(mag), MY(mag), MZ(mag)); 
     119    } 
     120  } 
     121  // Interpret polarization cross section. 
     122  double uu, dd, ud, du; 
     123  double cos_mspin, sin_mspin; 
     124  spins(up_frac_i, up_frac_f, &uu, &dd, &ud, &du); 
     125  SINCOS(-up_angle*M_PI_180, sin_mspin, cos_mspin); 
     126#endif 
    60127 
    61128  // Monodisperse computation 
    … …  
    74141    #endif 
    75142    for (int q_index=0; q_index < nq; q_index++) { 
     143#ifdef MAGNETIC 
     144      const double qx = q[2*q_index]; 
     145      const double qy = q[2*q_index+1]; 
     146      const double qsq = qx*qx + qy*qy; 
     147 
     148      // Constant across orientation, polydispersity for given qx, qy 
     149      double px, py, pz; 
     150      if (qsq > 1e-16) { 
     151        px = (qy*cos_mspin + qx*sin_mspin)/qsq; 
     152        py = (qy*sin_mspin - qx*cos_mspin)/qsq; 
     153        pz = 1.0; 
     154      } else { 
     155        px = py = pz = 0.0; 
     156      } 
     157 
     158      double scattering = 0.0; 
     159      if (uu > 1e-8) { 
     160        for (int mag=0; mag<NUM_MAGNETIC; mag++) { 
     161            const double perp = (qy*MX(mag) - qx*MY(mag)); 
     162            pvec[magnetic[mag]] = (values[magnetic[mag]+2] - perp*px)*uu; 
     163        } 
     164        scattering += CALL_IQ(q, q_index, local_values); 
     165      } 
     166      if (dd > 1e-8){ 
     167        for (int mag=0; mag<NUM_MAGNETIC; mag++) { 
     168            const double perp = (qy*MX(mag) - qx*MY(mag)); 
     169            pvec[magnetic[mag]] = (values[magnetic[mag]+2] + perp*px)*dd; 
     170        } 
     171        scattering += CALL_IQ(q, q_index, local_values); 
     172      } 
     173      if (ud > 1e-8){ 
     174        for (int mag=0; mag<NUM_MAGNETIC; mag++) { 
     175            const double perp = (qy*MX(mag) - qx*MY(mag)); 
     176            pvec[magnetic[mag]] = perp*py*ud; 
     177        } 
     178        scattering += CALL_IQ(q, q_index, local_values); 
     179        for (int mag=0; mag<NUM_MAGNETIC; mag++) { 
     180            pvec[magnetic[mag]] = MZ(mag)*pz*ud; 
     181        } 
     182        scattering += CALL_IQ(q, q_index, local_values); 
     183      } 
     184      if (du > 1e-8) { 
     185        for (int mag=0; mag<NUM_MAGNETIC; mag++) { 
     186            const double perp = (qy*MX(mag) - qx*MY(mag)); 
     187            pvec[magnetic[mag]] = perp*py*du; 
     188        } 
     189        scattering += CALL_IQ(q, q_index, local_values); 
     190        for (int mag=0; mag<NUM_MAGNETIC; mag++) { 
     191            pvec[magnetic[mag]] = -MZ(mag)*pz*du; 
     192        } 
     193        scattering += CALL_IQ(q, q_index, local_values); 
     194      } 
     195#else 
    76196      double scattering = CALL_IQ(q, q_index, local_values); 
     197#endif 
    77198      result[q_index] = (norm>0. ? scale*scattering/norm + background : background); 
    78199    } 
    … …  
    82203#if MAX_PD > 0 
    83204 
     205#if MAGNETIC 
     206  const double *pd_value = values+2+NPARS+3+3*NUM_MAGNETIC; 
     207#else 
    84208  const double *pd_value = values+2+NPARS; 
     209#endif 
    85210  const double *pd_weight = pd_value+details->pd_sum; 
    86211 
    … …  
    164289      #endif 
    165290      for (int q_index=0; q_index < nq; q_index++) { 
    166         const double scattering = CALL_IQ(q, q_index, local_values); 
     291#ifdef MAGNETIC 
     292        const double qx = q[2*q_index]; 
     293        const double qy = q[2*q_index+1]; 
     294        const double qsq = qx*qx + qy*qy; 
     295 
     296        // Constant across orientation, polydispersity for given qx, qy 
     297        double px, py, pz; 
     298        if (qsq > 1e-16) { 
     299          px = (qy*cos_mspin + qx*sin_mspin)/qsq; 
     300          py = (qy*sin_mspin - qx*cos_mspin)/qsq; 
     301          pz = 1.0; 
     302        } else { 
     303          px = py = pz = 0.0; 
     304        } 
     305 
     306        double scattering = 0.0; 
     307        if (uu > 1e-8) { 
     308          for (int mag=0; mag<NUM_MAGNETIC; mag++) { 
     309              const double perp = (qy*MX(mag) - qx*MY(mag)); 
     310              pvec[magnetic[mag]] = (values[magnetic[mag]+2] - perp*px)*uu; 
     311          } 
     312          scattering += CALL_IQ(q, q_index, local_values); 
     313        } 
     314        if (dd > 1e-8){ 
     315          for (int mag=0; mag<NUM_MAGNETIC; mag++) { 
     316              const double perp = (qy*MX(mag) - qx*MY(mag)); 
     317              pvec[magnetic[mag]] = (values[magnetic[mag]+2] + perp*px)*dd; 
     318          } 
     319          scattering += CALL_IQ(q, q_index, local_values); 
     320        } 
     321        if (ud > 1e-8){ 
     322          for (int mag=0; mag<NUM_MAGNETIC; mag++) { 
     323              const double perp = (qy*MX(mag) - qx*MY(mag)); 
     324              pvec[magnetic[mag]] = perp*py*ud; 
     325          } 
     326          scattering += CALL_IQ(q, q_index, local_values); 
     327          for (int mag=0; mag<NUM_MAGNETIC; mag++) { 
     328              pvec[magnetic[mag]] = MZ(mag)*pz*ud; 
     329          } 
     330          scattering += CALL_IQ(q, q_index, local_values); 
     331        } 
     332        if (du > 1e-8) { 
     333          for (int mag=0; mag<NUM_MAGNETIC; mag++) { 
     334              const double perp = (qy*MX(mag) - qx*MY(mag)); 
     335              pvec[magnetic[mag]] = perp*py*du; 
     336          } 
     337          scattering += CALL_IQ(q, q_index, local_values); 
     338          for (int mag=0; mag<NUM_MAGNETIC; mag++) { 
     339              pvec[magnetic[mag]] = -MZ(mag)*pz*du; 
     340          } 
     341          scattering += CALL_IQ(q, q_index, local_values); 
     342        } 
     343#else 
     344        double scattering = CALL_IQ(q, q_index, local_values); 
     345#endif 
    167346        result[q_index] += weight*scattering; 
    168347      } 
Note: See TracChangeset for help on using the changeset viewer.