Ticket #1243: six_shells_01.c

double lognormal(double r, double median, double sigma);
double Pq(double q, double t_c, double t_h, double t_p_in, double t_p_out, double rho_c, double rho_h, double rho_p_in, double rho_p_out, double rho_s_in, double rho_s_out, double R_current);
double Iq(double q, double t_c, double t_h, double t_p_in, double t_p_out, double rho_c, double rho_h, double rho_p_in, double rho_p_out, double rho_s_in, double rho_s_out, double R_median, double sigma_R, int n_pts, int n_sigs, int polyparam, double constraint);
double Iqxy(double qx, double qy, double t_c, double t_h, double t_p_in, double t_p_out, double rho_c, double rho_h, double rho_p_in, double rho_p_out, double rho_s_in, double rho_s_out, double R_median, double sigma_R, int n_pts, int n_sigs, int polyparam, double constraint);
double form_volume(double R_median, double sigma_R);

#define INVALID(v) (v.R_median<0, v.sigma_R<0, v.n_pts<0, v.n_pts>512, v.n_sigs<0, v.polyparam<1, v.polyparam>3, v.constraint<0)

#define max(a,b) \
   ({ __typeof__ (a) _a = (a); \
       __typeof__ (b) _b = (b); \
     _a > _b ? _a : _b; })

double lognormal(double r, double median, double sigma){ // this is only called for finite sigma, so no case sigma=0 covered
  const double sig = sigma/median;
  const double mu = log(median);
  if (r==0) {
    return 0;
  }
  else {
    return exp(-0.5*square((log(r)-mu)/sig))/(r*sig);
  }
}

double Pq(double q, double t_c, double t_h, double t_p_in, double t_p_out, double rho_c, double rho_h, double rho_p_in, double rho_p_out, double rho_s_in, double rho_s_out, double R_current){
  double R[8];
  double rho[8];
  
  R[3] = R_current - t_c/2; 
  R[4] = R_current + t_c/2;
  R[2] = R_current - t_c/2 - t_h;
  R[5] = R_current + t_c/2 + t_h;
  R[1] = R_current - t_c/2 - t_h - t_p_in;
  R[6] = R_current + t_c/2 + t_h + t_p_out;
  R[0] = 0.;
  
  R[7] = 0.;
  
  rho[0] = rho_s_in;
  rho[1] = rho_p_in;
  rho[2] = rho_h;
  rho[3] = rho_c;
  rho[4] = rho_h;
  rho[5] = rho_p_out;
  rho[6] = rho_s_out;
  
  rho[7] = 0;
  
  const double a0 = 2.8179403227e-6; // classical electron radius in nm
  
  double F = 0.;
  double V = 0.;
  for (int i=1;i<7;i++){
    V = M_4PI_3*cube(R[i]);
    F += V*(rho[i]-rho[i-1])*sas_3j1x_x(q*R[i]);
  }
  if (V==0){
    return 0;
  } else {
    return 1e6*a0*a0*F*F/V; // prefactor 1e6*a0^2 brings the intensity to mm^-1 if all radii are in nm and all densities are in e-/nm^3
  }
}

double Iq(double q, double t_c, double t_h, double t_p_in, double t_p_out, double rho_c, double rho_h, double rho_p_in, double rho_p_out, double rho_s_in, double rho_s_out, double R_median, double sigma_R, const int n_pts, const int n_sigs, int polyparam, double constraint){
  int i;
  double V;
  
  // size polydispersity distribution
  
  double r[512]; // this is an arbitrary number, one could use less
  double Pr[512]; // multiples of 2 probably give performance advantage
  
  const double Rmin = max(R_median-n_sigs*sigma_R,  t_p_in+t_h+t_c/2);
  const double Rmax = R_median+n_sigs*sigma_R;
  const double dR = (Rmax - Rmin)/((double)n_pts-1);
  
  double Iq = 0;
  if (sigma_R>0) {
    double norm = 0.;
    if(polyparam==1){ // conventional normalization
      for (i=0;i<n_pts;i++){
        r[i] = Rmin + i*dR;
        Pr[i] = lognormal(r[i],R_median,sigma_R);
        norm += Pr[i]*dR;
      }
      for (i=0;i<n_pts;i++){
        Pr[i] /= norm; // ignores constraint
      }
    } else if (polyparam==2){ // area normalization
      for (i=0;i<n_pts;i++){
        r[i] = Rmin + i*dR;
        Pr[i] = lognormal(r[i],R_median,sigma_R);
        norm += 4*M_PI*Pr[i]*(square(r[i]-t_c/2-t_h/2) + square(r[i]+t_c/2+t_h/2))*dR;
      }
      for (i=0;i<n_pts;i++){
        V = M_4PI_3*cube(r[i]);
        Pr[i] *= V*constraint/norm; // returns V*P(R)
      }
    } else if (polyparam==3){ // volume normalization
      for (i=0;i<n_pts;i++){
        r[i] = Rmin + i*dR;
        Pr[i] = lognormal(r[i],R_median,sigma_R);
        norm += M_4PI_3*Pr[i]*(cube(r[i]-t_c/2-t_h) - cube(r[i]-t_c/2-t_h-t_p_in) + cube(r[i]+t_c/2+t_h+t_p_out) - cube(r[i]+t_c/2+t_h))*dR;
      }
      for (i=0;i<n_pts;i++){
        V = M_4PI_3*cube(r[i]);
        Pr[i] *= V*constraint/norm; // returns V*P(R)
      }
    }
  // model function evaluation  
    for (i=0;i<n_pts;i++){
      Iq += Pr[i]*Pq(q, t_c, t_h, t_p_in, t_p_out, rho_c, rho_h, rho_p_in, rho_p_out, rho_s_in, rho_s_out, r[i])*dR;
    }    
  } else {
    Iq = Pq(q, t_c, t_h, t_p_in, t_p_out, rho_c, rho_h, rho_p_in, rho_p_out, rho_s_in, rho_s_out, R_median); // monodisperse case
  }
  return Iq;
}
double Iqxy(double qx, double qy, double t_c, double t_h, double t_p_in, double t_p_out, double rho_c, double rho_h, double rho_p_in, double rho_p_out, double rho_s_in, double rho_s_out, double R_median, double sigma_R, int n_pts, int n_sigs, int polyparam, double constraint){
  const double q = sqrt(qx*qx + qy*qy);
  return Iq(q, t_c, t_h, t_p_in, t_p_out, rho_c, rho_h, rho_p_in, rho_p_out, rho_s_in, rho_s_out, R_median, sigma_R, n_pts, n_sigs, polyparam, constraint);
}