yanera_reflectivity.c

Go to the documentation of this file.

/*==============================================================================
Copyright 2006 Thad Harroun, Kevin Yager

This file is part of "yanera 2.0".

"yanera 2.0" is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 3 of the License, or (at your option) any later version.

"yanera 2.0" is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.

You should have received a copy of the GNU General Public License along with "yanera 2.0"; if not, write to the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA 
==============================================================================*/ 
/*============================================================================*/
#define _CRT_SECURE_NO_DEPRECATE

#include "yanera.h"
#include "yanera_reflectivity.h"
#include "yanera_quadrature.h"
/*============================================================================*/
void yanera_write_reflectivity(yanera_container *yanera) 
{
  double            dq=0.0, q=0.0, r=0.0;
  int               i1, i2;
  FILE             *fp;
  char              filename[50];
  
  if (yanera->misc.fit_do == NO)
  {
    /*
     *  Find the step size in q.
     */
    if (yanera->misc.q_min == 0.0)  yanera->misc.q_min = 1.0E-5;
    dq = (yanera->misc.q_max - yanera->misc.q_min)/((double)yanera->misc.q_num);
  }
  
  for (i2=0; i2<yanera->number_of_models; i2++)
  {
    /* 
     * Make a filename for the reflectivity curve.
     */
    memset(filename, '\0', 50);
    if (yanera->data[i2].file_name == NULL) 
      { sprintf(filename, "result_%01d.ref", i2); }
    else sprintf(filename, "%s.ref", yanera->data[i2].file_name);
    
    /*
     * Break into slabs
     */
    if ((yanera->type == MODEL_COMPONENT) || \
        (yanera->type == MODEL_SLAB)      || \
        (yanera->type == MODEL_FUNCTION))
    { doQuadrature(&yanera->models_complete[i2], yanera); }
    
    /*
     * Open file and write to it.
     */
    fp = NULL;
    fp = fopen(filename,"w");
    if (fp == NULL) 
    {
      yanera_error("Can't open reflectivity file for writing", \
      __FILE__, __LINE__, yanera);
    }
    
    if (yanera->misc.fit_do == NO)
    {
      for (i1 = 0; i1 <= yanera->misc.q_num; i1++)
      {
        q = yanera->misc.q_min + dq * (double)i1;
        r = calcReflectivity(q, &yanera->data[i2].resolution, \
                                &yanera->models_complete[i2], yanera);
        fprintf(fp, "%g\t%g\n", q, r);
      }
    }
    else
    {
      for (i1 = 0; i1 < yanera->data[i2].n; i1++)
      {
        q = yanera->data[i2].q[i1];
        r = calcReflectivity(q, &yanera->data[i2].resolution, \
                                &yanera->models_complete[i2], yanera);
        fprintf(fp, "%g\t%g\n", q, r);
      }
    }
    
    fclose(fp);
  }
    
  return;
}
/*============================================================================*/
double calcReflectivity(double q, yanera_resolution *resolution, \
                        yanera_model *model, yanera_container *yanera)
{
  double tmp, value;

  tmp = yanera->parrattFunction(q, model, yanera);
 
  /* determine resolution for this q */
  value = 0.0;
  if ((resolution->type == RESOLUTION_CONSTANT) || \
      (resolution->type == RESOLUTION_RELATIVE))
  {
    value = resolution->value[0];
  }
  else if ((resolution->type == RESOLUTION_ARRAY_CONSTANT) || \
           (resolution->type == RESOLUTION_ARRAY_RELATIVE))
  {
    value = lookupResolution(q, resolution);
  }
  
  /* instrument resolution, skip if zero for speed */
  if (value != 0.0)
  {

    if ((resolution->type == RESOLUTION_CONSTANT) || \
        (resolution->type == RESOLUTION_ARRAY_CONSTANT))
    {
      tmp += 0.882*yanera->parrattFunction(q+(0.25*value), model, yanera);
      tmp += 0.882*yanera->parrattFunction(q-(0.25*value), model, yanera);
      tmp += 0.607*yanera->parrattFunction(q+(0.50*value), model, yanera);
      tmp += 0.607*yanera->parrattFunction(q-(0.50*value), model, yanera);
      tmp += 0.325*yanera->parrattFunction(q+(0.75*value), model, yanera);
      tmp += 0.325*yanera->parrattFunction(q-(0.75*value), model, yanera);
      tmp += 0.135*yanera->parrattFunction(q+(1.00*value), model, yanera);
      tmp += 0.135*yanera->parrattFunction(q-(1.00*value), model, yanera);
    }
    else if ((resolution->type == RESOLUTION_RELATIVE) || \
             (resolution->type == RESOLUTION_ARRAY_RELATIVE))
    {
      tmp += 0.882*yanera->parrattFunction(q*(1.0+(0.25*value)), model, yanera);
      tmp += 0.882*yanera->parrattFunction(q*(1.0-(0.25*value)), model, yanera);
      tmp += 0.607*yanera->parrattFunction(q*(1.0+(0.50*value)), model, yanera);
      tmp += 0.607*yanera->parrattFunction(q*(1.0-(0.50*value)), model, yanera);
      tmp += 0.325*yanera->parrattFunction(q*(1.0+(0.75*value)), model, yanera);
      tmp += 0.325*yanera->parrattFunction(q*(1.0-(0.75*value)), model, yanera);
      tmp += 0.135*yanera->parrattFunction(q*(1.0+(1.00*value)), model, yanera);
      tmp += 0.135*yanera->parrattFunction(q*(1.0-(1.00*value)), model, yanera);
    }
    tmp /= (4.898);
  }
  
  /* background */
  if (model->background > -1)
    tmp += (yanera->parameters.p[model->background] * 1.0E-06);
    
 return tmp;
  
}
/*============================================================================*/
double calcParrattFromSlabs(double q, \
                            yanera_model *model, \
                            yanera_container *yanera)
{
  short        i;
  double       kv, tmpr, tmpi, result;
  gsl_complex  kzn, kznp1;
  gsl_complex  fr, R;
  gsl_complex  tmp1, tmp2;
  gsl_complex  argm, numtr, denom;
  
  if (q <= 0.0) return 1.0;
  
  /*
   * k_v: wave vector normal to the interface in vacuum.
   * k_0 = n_0*k_v \approx q/2. 
   */
  tmpr = (q * q / 4.0) + (4.0 * M_PI * model->slabs.rsld[0] * 1E-6);
  tmpi = (4.0 * M_PI * model->slabs.isld[0] * 1E-6);
  GSL_SET_COMPLEX(&tmp1, tmpr, tmpi);
  kv  = gsl_complex_abs(gsl_complex_sqrt(tmp1));
  
  /*
   * Zero the refecltivty variables: 
   * R = total reflectivity
   * fr = Fresnel reflectivity from any interface.
   */ 
  GSL_SET_COMPLEX(&R, 0.0, 0.0);
  GSL_SET_COMPLEX(&fr, 0.0, 0.0);
  
  /*
   * FORMULA INDEXING OF THE SLABS
   * If there are N slabs, infinite fronting material is slab 0, 
   * infinite backing material is slab N+1.
   * First reflection, from interface between slab N and N+1(bulk),
   * is Fresnel reflection only. 
   *
   * C INDEXING OF THE SLABS
   * The C index of n slabs (counting fronting and backing material)
   * is 0->n-1, so this first interface is n-2 and n-1.
   *
   * kzn   : wave vector in slab n
   * kznp1 : wave vector in slab n+1
   */
  i = model->slabs.number_of_slabs-2;
  {
    tmpr = (kv * kv)-(4.0 * M_PI * model->slabs.rsld[i] * 1E-6);
    tmpi = (4.0 * M_PI * model->slabs.isld[i] * 1E-6);
    GSL_SET_COMPLEX(&tmp1, tmpr, tmpi);
    kzn   = gsl_complex_sqrt(tmp1);
    
    tmpr = (kv * kv)-(4.0 * M_PI * model->slabs.rsld[i+1] * 1E-6);
    tmpi = (4.0 * M_PI * model->slabs.isld[i+1] * 1E-6);
    GSL_SET_COMPLEX(&tmp1, tmpr, tmpi);
    kznp1   = gsl_complex_sqrt(tmp1);
    
    /* Fresnel */
    tmp1 = gsl_complex_sub(kzn, kznp1);
    tmp2 = gsl_complex_add(kzn, kznp1);
    R    = gsl_complex_div(tmp1, tmp2);
  }
  
  /* 
   *  Subsequent reflections: End on the reflection from slab 0->1.
   */
  for (i=model->slabs.number_of_slabs-3; i>=0; i--)
  {
    GSL_SET_COMPLEX(&kznp1, GSL_REAL(kzn), GSL_IMAG(kzn));
  
    tmpr = (kv * kv)-(4.0 * M_PI * model->slabs.rsld[i] * 1E-6);
    tmpi = (4.0 * M_PI * model->slabs.isld[i] * 1E-6);
    GSL_SET_COMPLEX(&tmp1, tmpr, tmpi);
    kzn   = gsl_complex_sqrt(tmp1);
    
    /* 
     * Save one sqrt(), use the fact that this already calculated during the 
     * previous iteration!
    tmpr = (kv * kv)-(4.0 * M_PI * model->slabs.rsld[i+1] * 1E-6);
    tmpi = (4.0 * M_PI * model->slabs.isld[i+1] * 1E-6);
    GSL_SET_COMPLEX(&tmp1, tmpr, tmpi);
    kznp1   = gsl_complex_sqrt(tmp1);
    */
    
    /* Fresnel */
    tmp1 = gsl_complex_sub(kzn, kznp1);
    tmp2 = gsl_complex_add(kzn, kznp1);
    fr   = gsl_complex_div(tmp1, tmp2);
  
    /* argm = R_n+1 * exp(2 i k_n+1 d_n+1) */
    tmp1 = gsl_complex_mul_imag(kznp1, 2.0*model->slabs.thik[i+1]);
    tmp2 = gsl_complex_exp(tmp1);
    argm = gsl_complex_mul(R, tmp2);

    /* numtr = fr_n,n+1 + R_n+1 * exp(2 i k_n+1 d_n+1) */
    numtr = gsl_complex_add(fr, argm);
        
    /* denom = 1 + fr_n,n+1 * R_n+1 * exp(2 i k_n+1 d_n+1) */
    denom = gsl_complex_add_real(gsl_complex_mul(fr, argm), 1.0);
    
    R = gsl_complex_div(numtr, denom);
  }
  
  /* Grand result: |R|^2 */
  result = gsl_complex_abs2(R);
  
  return result;
}
/*============================================================================*/
double calcParrattFromModel(double q, \
                            yanera_model *model, \
                            yanera_container *yanera)
{
  double        kv, tmpr, tmpi, result, rsld;
  gsl_complex   kzn, kznp1;
  gsl_complex   fr, r, R;
  gsl_complex   tmp1, tmp2, tmp3;
  gsl_complex   argm, numtr, denom;
  yanera_layer *layer;
  
  if (q < 0) return 1.0;
  
  /*
   * k_v: wave vector normal to the interface in vacuum.
   * k_0 = n_0*k_v \approx q/2. 
   */
  layer = model->first_layer;
  rsld = yanera->parameters.p[layer->rsld];
  if      (layer->polarized == POLARIZED_UP) 
    rsld += yanera->parameters.p[layer->rmag];
  else if (layer->polarized == POLARIZED_DOWN) 
    rsld -= yanera->parameters.p[layer->rmag];
  tmpr = (q * q / 4.0) + (4.0 * M_PI * rsld * 1E-6);
  if (layer->isld > -1)  
  { tmpi = (4.0 * M_PI * yanera->parameters.p[layer->isld] * 1E-6); }
  else tmpi = 0.0;
  GSL_SET_COMPLEX(&tmp1, tmpr, tmpi);
  kv  = gsl_complex_abs(gsl_complex_sqrt(tmp1));
  
  /*
   * Zero the refecltivty variables: 
   * R = total reflectivity
   * fr = Fresnel reflectivity from any interface.
   */ 
  GSL_SET_COMPLEX(&R, 0.0, 0.0);
  GSL_SET_COMPLEX(&fr, 0.0, 0.0);
  GSL_SET_COMPLEX(&r, 0.0, 0.0);
  
  /*
   * FORMULA INDEXING OF THE SLABS
   * If there are N layers, infinite fronting material is layer 0, 
   * infinite backing material is layer N+1.
   * First reflection, from interface between layers N and N+1(bulk),
   * is Fresnel reflection only. 
   *
   * C INDEXING OF THE SLABS
   * The C index of n layers (counting fronting and backing material)
   * is 0->n-1, so this interface is n-2 and n-1.
   *
   * kzn   : wave vector in slab n
   * kznp1 : wave vector in slab n+1
   */
  layer = model->last_layer->prev;
  {
    rsld = yanera->parameters.p[layer->rsld];
    if      (layer->polarized == POLARIZED_UP) 
      rsld += yanera->parameters.p[layer->rmag];
    else if (layer->polarized == POLARIZED_DOWN) 
      rsld -= yanera->parameters.p[layer->rmag];
    tmpr = (kv * kv)-(4.0 * M_PI * rsld * 1E-6);
    if (layer->isld > -1)  
    { tmpi = (4.0 * M_PI * yanera->parameters.p[layer->isld] * 1E-6); }
    else tmpi = 0.0;
    GSL_SET_COMPLEX(&tmp1, tmpr, tmpi);
    kzn   = gsl_complex_sqrt(tmp1);
    
    rsld = yanera->parameters.p[layer->next->rsld];
    if      (layer->next->polarized == POLARIZED_UP) 
      rsld += yanera->parameters.p[layer->next->rmag];
    else if (layer->next->polarized == POLARIZED_DOWN) 
      rsld -= yanera->parameters.p[layer->next->rmag];
    tmpr = (kv * kv)-(4.0 * M_PI * rsld * 1E-6);
    if (layer->next->isld > -1)  
    { tmpi = (4.0 * M_PI * yanera->parameters.p[layer->next->isld] * 1E-6); }
    else tmpi = 0.0;
    GSL_SET_COMPLEX(&tmp1, tmpr, tmpi);
    kznp1   = gsl_complex_sqrt(tmp1);
    
    /* Fresnel */
    tmp1 = gsl_complex_sub(kzn, kznp1);
    tmp2 = gsl_complex_add(kzn, kznp1);
    fr    = gsl_complex_div(tmp1, tmp2);
    
    /* r = Fresnel * roughness */
    tmp1 = gsl_complex_mul(kzn,kznp1);
    tmpr = -2.0 * yanera->parameters.p[layer->next->sigz] * \
                  yanera->parameters.p[layer->next->sigz];
    tmp2 = gsl_complex_mul_real(tmp1, tmpr);
    tmp3 = gsl_complex_exp(tmp2);
    R    = gsl_complex_mul(fr, tmp3);
  }
  
  /* 
   *  Subsequent reflections: End on the reflection from slab 0->1.
   */
  do
  {
    layer = layer->prev;
    
    GSL_SET_COMPLEX(&kznp1, GSL_REAL(kzn), GSL_IMAG(kzn));
    
    rsld = yanera->parameters.p[layer->rsld];
    if      (layer->polarized == POLARIZED_UP) 
      rsld += yanera->parameters.p[layer->rmag];
    else if (layer->polarized == POLARIZED_DOWN) 
      rsld -= yanera->parameters.p[layer->rmag];
    tmpr = (kv * kv)-(4.0 * M_PI * rsld * 1E-6);
    if (layer->isld > -1)  
    { tmpi = (4.0 * M_PI * yanera->parameters.p[layer->isld] * 1E-6); }
    else tmpi = 0.0;
    GSL_SET_COMPLEX(&tmp1, tmpr, tmpi);
    kzn   = gsl_complex_sqrt(tmp1);
    
    /* 
     * Save one sqrt(), use the fact that this already calculated during the 
     * previous iteration!
    rsld = yanera->parameters.p[layer->next->rsld];
    if      (layer->next->polarized == POLARIZED_UP) 
      rsld += yanera->parameters.p[layer->next->rmag];
    else if (layer->next->polarized == POLARIZED_DOWN) 
      rsld -= yanera->parameters.p[layer->next->rmag];
    tmpr = (kv * kv)-(4.0 * M_PI * rsld * 1E-6);
    if (layer->next->isld > -1)  
    { tmpi = (4.0 * M_PI * yanera->parameters.p[layer->next->isld] * 1E-6); }
    else tmpi = 0.0;
    GSL_SET_COMPLEX(&tmp1, tmpr, tmpi);
    kznp1   = gsl_complex_sqrt(tmp1);
    */
    
    /* Fresnel */
    tmp1 = gsl_complex_sub(kzn, kznp1);
    tmp2 = gsl_complex_add(kzn, kznp1);
    fr   = gsl_complex_div(tmp1, tmp2);
    
    /* r = Fresnel * roughness */
    tmp1 = gsl_complex_mul(kzn,kznp1);
    tmpr = -2.0 * yanera->parameters.p[layer->next->sigz] * \
                  yanera->parameters.p[layer->next->sigz];
    tmp2 = gsl_complex_mul_real(tmp1, tmpr);
    tmp3 = gsl_complex_exp(tmp2);
    r    = gsl_complex_mul(fr, tmp3);
  
    /* argm = R_n+1 * exp(2 i k_n+1 d_n+1) */
    tmp1 = gsl_complex_mul_imag(kznp1, \
                          2.0*yanera->parameters.p[layer->next->thik]);
    tmp2 = gsl_complex_exp(tmp1);
    argm = gsl_complex_mul(R, tmp2);

    /* numtr = r_n,n+1 + R_n+1 * exp(2 i k_n+1 d_n+1) */
    numtr = gsl_complex_add(r, argm);
        
    /* denom = 1 + r_n,n+1 * R_n+1 * exp(2 i k_n+1 d_n+1) */
    denom = gsl_complex_add_real(gsl_complex_mul(r, argm), 1.0);
    
    R = gsl_complex_div(numtr, denom);
    
  } while (layer != model->first_layer);
  
  /* Grand result: |R|^2 */
  result = gsl_complex_abs2(R);
  
  return result;
}
/*============================================================================*/
double lookupResolution(double  q, yanera_resolution *resolution)
{
  int i=0;
  /* 
   *  Walk through the array of q-values in reverse until you walk past
   *  the target value 
   */
  for (i = resolution->n-1; i>=0; i--)
  {
    if (resolution->q[i] <= q) return resolution->value[i];
  }
  
  /* If no resolution is specified for low-q data, use 0 */
  return 0.0;
 }

---