yanera_quadrature.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_quadrature.h"
#include "yanera_profile.h"
/* -------------------------------------------------------------------------- */
double slabs_trap(yanera_container *p, yanera_model *m, unsigned short i)
{
  double tmp;
  
  tmp = (m->slabs.zpos[i+1] - m->slabs.zpos[i]) * \
        (p->profileFunctionReal(m->slabs.zpos[i+1], m, p)+ \
         p->profileFunctionReal(m->slabs.zpos[i],   m, p))/2.0;
  return tmp;
}
/* -------------------------------------------------------------------------- */
double slabs_midp(yanera_container *p, yanera_model *m, unsigned short i)
{
  double tmp1, tmp2;
  tmp1 = (m->slabs.zpos[i+1] + m->slabs.zpos[i])/2.0;
  tmp2 = (m->slabs.zpos[i+1] - m->slabs.zpos[i]) * \
          p->profileFunctionReal(tmp1, m, p);
  return tmp2;
}
/* -------------------------------------------------------------------------- */
void slabs_insert(yanera_model *model, unsigned short j, double np)
{
  unsigned short i;
  
  if ((model->slabs.number_of_slab_edges+1) >= USHRT_MAX) 
    yanera_error("Too many slabs.", __FILE__, __LINE__, NULL);
  
  for (i=model->slabs.number_of_slab_edges-1; i>j; i--)
  {
    model->slabs.zpos[i+1] = model->slabs.zpos[i];
  }
  model->slabs.zpos[j+1]              = np;
  model->slabs.number_of_slab_edges  += 1;
  model->slabs.number_of_slabs       += 1;
    
  return;
}
/* -------------------------------------------------------------------------- */
void slabs_zero(yanera_model *model, double start, double end)
{
  model->slabs.number_of_slab_edges = 1;
  model->slabs.number_of_slabs      = 2;
  model->slabs.zpos[0]              = start;
  model->slabs.thik[0]              = 0.0;
  
  slabs_insert(model, 0, end);
  
  return;
}
/* -------------------------------------------------------------------------- */
void slabs_doit(yanera_model *m, yanera_container *p)
{
  unsigned short i;
  double         tmp, zpos;
  
  i = 0;
  while (i != m->slabs.number_of_slab_edges-1)
  {
    tmp = fabs(slabs_trap(p,m,i)-slabs_midp(p,m,i));
    if (tmp > p->misc.quadrature_error)
    {
      zpos = (m->slabs.zpos[i] + m->slabs.zpos[i+1]) / 2.0;
      slabs_insert(m, i, zpos);
    }
    else
    {
      /* 
       * This bit here is to re-start the adaptive quadrature in the event
       * it makes a slab that is too big.
       * This can happen in periodic structures, when the quadarutre stretches
       * across an number of layers that happen to have the same rsld. In that
       * case, the integral difference (tmp) is very small, and the
       * subdivision stops.
       * 
       * First, if we reached here, the intregral difference (tmp) is smaller
       * than the tolerance. Normally, we would move to the next interval, and
       * consider the subdivision of this interval small enough.
       * But if the interval is bigger than quadrature_thik, 
       * then it really isn't small enough: divide it again 
       * and continue dividing as needed.
       * 
       * This puts an inherent limit on the total film thickness to about
       * USHRT_MAX * quadrature_thik. Is this enough? 
       */
      zpos = (m->slabs.zpos[i+1] - m->slabs.zpos[i]);
      if (zpos > p->misc.quadrature_thik)
      {
        zpos = m->slabs.zpos[i] + zpos/2.0;
        slabs_insert(m, i, zpos);
      }
      else
      {
        i += 1;
      }
    }
  }
 
  return;
}
/* -------------------------------------------------------------------------- */
void doQuadrature(yanera_model *model, yanera_container *yanera)
{
  double         start, end, zpos;
  short          i;
  
  /* Get the min/max of the rsld profile */
  minMaxProfile(&start, &end, model, yanera);
    
  /*
   *  Do the quadrature
   *  After, we have the edges of all the slabs.
   */
  slabs_zero(model, start, end);
  slabs_doit(model, yanera);
    
  /*
   * Since we store the edge and thickness together, define the edge (z)
   * of the slab as the high-z side.
   * 
   * Index slabs as : 0 to slabs.number_of_slabs-1
   * Index z points found by the quadrature (the slab edges) as:
   * 0 to slabs.number_of_slab_edges-1.
   *
   * Let : slabs.number_of_slabs = number_of_slab_edges + 1
   *
   * For slab 0 and n-1 there is no thickness.
   * For slab 0, z is the 'start' (see above) of the profile.
   * For slab n-2, z value is the 'end' (see above) of the profile.
   * For slab n-1, there is no z.
   */
  for (i=model->slabs.number_of_slabs-2; i>0; i--)
  {
    model->slabs.thik[i] = model->slabs.zpos[i] - model->slabs.zpos[i-1];
                                     
    zpos = (model->slabs.zpos[i] + model->slabs.zpos[i-1]) / 2.0;
    model->slabs.rsld[i] = yanera->profileFunctionReal(zpos, model, yanera);
    model->slabs.isld[i] = yanera->profileFunctionImag(zpos, model, yanera);
  }
  /* Set the rsld of (infinite fronting material) slab 0 */
  model->slabs.rsld[0] = yanera->profileFunctionReal(start, model, yanera);
  model->slabs.isld[0] = yanera->profileFunctionImag(start, model, yanera);
  /* Set the rsld of (infinite backing material) slab n-1 */
  model->slabs.rsld[model->slabs.number_of_slabs-1] = \
                   yanera->profileFunctionReal(end, model, yanera);
  model->slabs.isld[model->slabs.number_of_slabs-1] = \
                   yanera->profileFunctionImag(end, model, yanera);
   
  return;
}
/*============================================================================*/
void yanera_write_slabs(yanera_container *yanera) 
{
  double           zpos;
  char             filename[50];
  short            i2, i1;
  yanera_model    *model;
  FILE            *fp;
  
  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.slb", i2); }
    else sprintf(filename, "%s.slb", yanera->data[i2].file_name);

    /*
     * Break into slabs
     */
    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);
    }
    /*
     * Use the pointer "model" to keep the code simple:
     * else &yanera->models_complete[i2] would be appended to everything below.
     */
    model = &yanera->models_complete[i2];
    /*
     * Write the first slab, which is the infinite bakcing, as extending
     * one angstrom to higher z.
     */
    zpos = model->slabs.zpos[model->slabs.number_of_slabs-2];
    fprintf(fp, "%f\t%f\n", zpos+1.0, \
          model->slabs.rsld[model->slabs.number_of_slabs-1]);
    fprintf(fp, "%f\t%f\n", zpos, \
          model->slabs.rsld[model->slabs.number_of_slabs-1]);
    /*
     * Write the middle slabs. Start at the point last left off.
     */
    for (i1=model->slabs.number_of_slabs-2; i1>0; i1--)
    {
      fprintf(fp, "%f\t%f\n", zpos, model->slabs.rsld[i1]);
      zpos -= model->slabs.thik[i1];
      fprintf(fp, "%f\t%f\n", zpos, model->slabs.rsld[i1]);
    }
    /*
     * Write the final, infinite fronting slab as one angstrom thick.
     */
    fprintf(fp, "%f\t%f\n", zpos,     model->slabs.rsld[0]);    
    fprintf(fp, "%f\t%f\n", zpos-1.0, model->slabs.rsld[0]);
    /* --- */ 
    fclose(fp);
  }
   
  return;    
}       

---