yanera_profile.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_profile.h"
#include "yanera_postfix.h"
/*============================================================================*/
void yanera_write_profile(yanera_container *yanera) 
{
  double           minz, maxz, dz, z, rsld, isld;
  short            i1, i2;
  FILE            *fp;
  char             filename[50];
  
  /*
   * Loop over the models. Create a filename and call the
   * function that writes the profile to file.
   */
  for (i2=0; i2<yanera->number_of_models; i2++)
  {
    /*
     * Set the filename of the profile.
     */
    memset(filename, '\0', 50);
    if (yanera->data[i2].file_name == NULL) 
    { sprintf(filename, "result_%01d.pro", i2); }
    else sprintf(filename, "%s.pro", yanera->data[i2].file_name);

    /*
     * Get the min and max z; calculate the step size
     */
    minMaxProfile(&minz, &maxz, &yanera->models_complete[i2], yanera);
    dz = (maxz - minz)/(double)yanera->misc.z_num;
  
    /*
     * Open the file for writing
     */
    fp = NULL;
    fp = fopen(filename,"w");
    if (fp == NULL) 
    {
      yanera_error("Can't open profile file for writing.", \
                   __FILE__, __LINE__, yanera);
    }
    /*
     * Write to file and close
     */
    for (i1 = 0; i1 <= yanera->misc.z_num; i1++)
    {
      z = minz + dz * (double)i1;
      rsld = yanera->profileFunctionReal(z, &yanera->models_complete[i2], \
                                         yanera);
      isld = yanera->profileFunctionImag(z, &yanera->models_complete[i2], \
                                         yanera);
      fprintf(fp, "%f\t%f\t%f\n", z, rsld, isld);
    }
    fclose(fp);
  }
    
  return;
}
/*============================================================================*/
double layerProfileReal(double z, yanera_model *model, yanera_container *yanera)
{
  double        thik, sld, tmp1, rho1, rho2;
  yanera_layer *layer1, *layer2;

  /* ------------------------------------
   * The first layer is the layer the neutrons enter from.
   */
  layer1 = model->first_layer;
  sld = yanera->parameters.p[layer1->rsld];
  if (layer1->polarized == POLARIZED_UP) 
    sld += yanera->parameters.p[layer1->rmag];
  else if (layer1->polarized == POLARIZED_DOWN) 
    sld -= yanera->parameters.p[layer1->rmag];
    
    
  /* ------------------------------------
   * Loop over the layers.
   */
  layer2 = layer1;
  layer1 = layer1->next;
  thik = 0.0;
  /* ---- */
  while (layer1 != NULL)
  {
    
    rho1   = yanera->parameters.p[layer1->rsld];
    if (layer1->polarized == POLARIZED_UP) 
       rho1 += yanera->parameters.p[layer1->rmag];
    else if (layer1->polarized == POLARIZED_DOWN) 
       rho1 -= yanera->parameters.p[layer1->rmag];
   
    
    rho2   = yanera->parameters.p[layer2->rsld];
    if (layer2->polarized == POLARIZED_UP) 
      rho2 += yanera->parameters.p[layer2->rmag];
    else if (layer2->polarized == POLARIZED_DOWN) 
      rho2 -= yanera->parameters.p[layer2->rmag];  

    if (yanera->parameters.p[layer1->sigz] == 0)
    {
      if (z >= thik) sld += (rho1 - rho2);
    }
    else
    {
      tmp1 = (z-thik)/(M_SQRT2*fabs(yanera->parameters.p[layer1->sigz]));
      sld += (rho1 - rho2) * (1.0 + gsl_sf_erf(tmp1)) / 2.0;
    }

    /* ------------------------------------ 
     * This *if* is because the last layer does not have a thik parameter
     * and without it, we would be accessing the parameter array 
     * outside its bounds.
     */
    if (layer1->thik > -1)
      thik += yanera->parameters.p[layer1->thik];
    
    layer1 = layer1->next;
    layer2 = layer2->next;
  }

  return sld;
}
/*============================================================================*/
double componentProfileReal(double z, yanera_model *model, \
                            yanera_container *yanera)
{
  double        zpos, sld, tmp1, tmp2, tmpl, tmph, tmpsld;
  yanera_layer *layer;

  sld = 0.0;
  /*=========================
   * The first layer is the layer the neutrons enter from. 
   * The simple form is a backwards facing error function
   * rho*(1-erf(z/(sqrt(2)*sigma)))
   */
  layer = model->first_layer;
  tmpsld = yanera->parameters.p[layer->rsld];
  if (layer->polarized == POLARIZED_UP) 
    tmpsld += yanera->parameters.p[layer->rmag];
  else if (layer->polarized == POLARIZED_DOWN) 
    tmpsld -= yanera->parameters.p[layer->rmag];
  zpos = 0.0;
  if (yanera->parameters.p[layer->sigz] == 0.0) 
  {
    if (z < zpos) sld += tmpsld;
  }
  else
  {
    tmp1 = (z-zpos)/(M_SQRT2*fabs(yanera->parameters.p[layer->sigz]));
    sld  += tmpsld * (1.0 - gsl_sf_erf(tmp1)) / 2.0;
  } 
    
  /*=========================
   * Loop over the next layers.
   */
  layer=layer->next; 
  while (layer->next != NULL)
  {
    tmpsld = yanera->parameters.p[layer->rsld];
    if (layer->polarized == POLARIZED_UP) 
      tmpsld += yanera->parameters.p[layer->rmag];
    else if (layer->polarized == POLARIZED_DOWN) 
      tmpsld -= yanera->parameters.p[layer->rmag];
    zpos = yanera->parameters.p[layer->cntr];

    if (layer->type == LAYER_BOX)
    {
      tmpl = zpos - yanera->parameters.p[layer->thik]/2.0; 
      tmph = zpos + yanera->parameters.p[layer->thik]/2.0; 
      /* zero roughness : erf becomes box function */
      if (yanera->parameters.p[layer->sigz]== 0.0)
      {
        if ((z>tmpl)&&(z<tmph)) sld += tmpsld;
      }
      else
      {
        tmp1 = (z-tmpl)/(M_SQRT2*fabs(yanera->parameters.p[layer->sigz]));
        tmp2 = (z-tmph)/(M_SQRT2*fabs(yanera->parameters.p[layer->sigz]));
        sld +=  tmpsld * \
                (gsl_sf_erf(tmp1) - gsl_sf_erf(tmp2)) / 2.0;
      }
      
    }
    else if (layer->type == LAYER_GAUSSIAN)
    {
      tmp1 = (zpos - z) / yanera->parameters.p[layer->sigz];
      sld += tmpsld * exp(-1.0 * tmp1 * tmp1);
    }
    else 
    {
      yanera_error("Unknown component type.", __FILE__, __LINE__, yanera);
    }
    layer = layer->next;
  }
  
  /*=========================
   * The last layer 
   * The simple form is a forwards facing error function
   * rho*(1+erf((z-cntr)/(sqrt(2)*sigma)))
   */
  tmpsld = yanera->parameters.p[layer->rsld];
  if (layer->polarized == POLARIZED_UP) 
    tmpsld += yanera->parameters.p[layer->rmag];
  else if (layer->polarized == POLARIZED_DOWN) 
    tmpsld -= yanera->parameters.p[layer->rmag];
  zpos = yanera->parameters.p[layer->cntr];
  if (yanera->parameters.p[layer->sigz] == 0) 
  {
    if (z >= zpos) sld += tmpsld;
  }
  else
  {
    tmp1 = (z-zpos)/(M_SQRT2*fabs(yanera->parameters.p[layer->sigz]));
    sld  += tmpsld * (1.0 + gsl_sf_erf(tmp1)) / 2.0;
  }    
  return sld;
}
/*============================================================================*/
double functionProfileReal(double z, yanera_model *model, \
                           yanera_container *yanera)
{
  double        thik, sld, tmp1, rho1, rho2;
  yanera_layer *layer1, *layer2;

  /* ------------------------------------
   * The first layer is the layer the neutrons enter from.
   */
  layer1 = model->first_layer;
  sld = yanera->parameters.p[layer1->rsld];
  if (layer1->polarized == POLARIZED_UP) 
    sld += yanera->parameters.p[layer1->rmag];
  else if (layer1->polarized == POLARIZED_DOWN) 
    sld -= yanera->parameters.p[layer1->rmag];
    
    
  /* ------------------------------------
   * Loop over the layers.
   */
  layer2 = layer1;
  layer1 = layer1->next;
  thik = 0.0;
  /* ---- */
  while (layer1 != NULL)
  {
    rho1 = 0;
    if (layer1->type == LAYER_FUNCTION)
    {
      rho1 = evaluate_postfix(layer1->func_postfix, z, yanera->parameters.p);
    }
    else
    {
      rho1   = yanera->parameters.p[layer1->rsld];
      if (layer1->polarized == POLARIZED_UP) 
        rho1 += yanera->parameters.p[layer1->rmag];
      else if (layer1->polarized == POLARIZED_DOWN) 
        rho1 -= yanera->parameters.p[layer1->rmag];
    }
    
    rho2 = 0;
    if (layer2->type == LAYER_FUNCTION)
    {
      rho2 = evaluate_postfix(layer2->func_postfix, z, yanera->parameters.p);
    }
    else
    {
      rho2   = yanera->parameters.p[layer2->rsld];
      if (layer2->polarized == POLARIZED_UP) 
        rho2 += yanera->parameters.p[layer2->rmag];
      else if (layer2->polarized == POLARIZED_DOWN) 
        rho2 -= yanera->parameters.p[layer2->rmag];
    }    

    if (yanera->parameters.p[layer1->sigz] == 0)
    {
      if (z >= thik) sld += (rho1 - rho2);
    }
    else
    {
      tmp1 = (z-thik)/(M_SQRT2*fabs(yanera->parameters.p[layer1->sigz]));
      sld += (rho1 - rho2) * (1.0 + gsl_sf_erf(tmp1)) / 2.0;
    }

    /* ------------------------------------ 
     * This *if* is because the last layer does not have a thik parameter
     * and without it, we would be accessing the parameter array 
     * outside its bounds.
     */
    if (layer1->thik > -1)
      thik += yanera->parameters.p[layer1->thik];
    
    layer1 = layer1->next;
    layer2 = layer2->next;
  }
  
  return sld;
}
/*============================================================================*/
double layerProfileImag(double z, yanera_model *model, yanera_container *yanera)
{
  double        thik, sld, tmp1, rho1, rho2;
  yanera_layer *layer1, *layer2;

  /* ------------------------------------
   * The first layer is the layer the neutrons enter from.
   */
  layer1 = model->first_layer;
  if (layer1->isld > -1) sld = yanera->parameters.p[layer1->isld];
  else sld = 0.0;
  
  /* ------------------------------------
   * Loop over the layers.
   */
  layer2 = layer1;
  layer1 = layer1->next;
  thik = 0.0;
  /* ---- */
  while (layer1 != NULL)
  {
    if (layer1->isld > -1) rho1 = yanera->parameters.p[layer1->isld];
    else rho1 = 0.0;
    
    if (layer2->isld > -1) rho2   = yanera->parameters.p[layer2->isld];
    else rho2 = 0.0;  

    if (yanera->parameters.p[layer1->sigz] == 0)
    {
      if (z >= thik) sld += (rho1 - rho2);
    }
    else
    {
      tmp1 = (z-thik)/(M_SQRT2*fabs(yanera->parameters.p[layer1->sigz]));
      sld += (rho1 - rho2) * (1.0 + gsl_sf_erf(tmp1)) / 2.0;
    }

    /* ------------------------------------ 
     * This *if* is because the last layer does not have a thik parameter
     * and without it, we would be accessing the parameter array 
     * outside its bounds.
     */
    if (layer1->thik > -1)
      thik += yanera->parameters.p[layer1->thik];
    
    layer1 = layer1->next;
    layer2 = layer2->next;
  }

  return sld;
}
/*============================================================================*/
double componentProfileImag(double z, yanera_model *model, \
                            yanera_container *yanera)
{
  double        zpos, sld, tmp1, tmp2, tmpl, tmph, tmpsld;
  yanera_layer *layer;

  sld = 0.0;
  /*=========================
   * The first layer is the layer the neutrons enter from. 
   * The simple form is a backwards facing error function
   * rho*(1-erf(z/(sqrt(2)*sigma)))
   */
  layer = model->first_layer;
  if (layer->isld > -1) tmpsld = yanera->parameters.p[layer->isld];
  else tmpsld = 0.0;
  zpos = 0.0;
  if (yanera->parameters.p[layer->sigz] == 0.0) 
  {
    if (z < zpos) sld += tmpsld;
  }
  else
  {
    tmp1 = (z-zpos)/(M_SQRT2*fabs(yanera->parameters.p[layer->sigz]));
    sld  += tmpsld * (1.0 - gsl_sf_erf(tmp1)) / 2.0;
  } 
    
  /*=========================
   * Loop over the next layers.
   */
  layer=layer->next; 
  while (layer->next != NULL)
  {
    if (layer->isld > -1) tmpsld = yanera->parameters.p[layer->isld];
    else tmpsld = 0.0;
    zpos = yanera->parameters.p[layer->cntr];

    if (layer->type == LAYER_BOX)
    {
      tmpl = zpos - yanera->parameters.p[layer->thik]/2.0; 
      tmph = zpos + yanera->parameters.p[layer->thik]/2.0; 
      /* zero roughness : erf becomes box function */
      if (yanera->parameters.p[layer->sigz]== 0.0)
      {
        if ((z>tmpl)&&(z<tmph)) sld += tmpsld;
      }
      else
      {
        tmp1 = (z-tmpl)/(M_SQRT2*fabs(yanera->parameters.p[layer->sigz]));
        tmp2 = (z-tmph)/(M_SQRT2*fabs(yanera->parameters.p[layer->sigz]));
        sld +=  tmpsld * \
                (gsl_sf_erf(tmp1) - gsl_sf_erf(tmp2)) / 2.0;
      }
      
    }
    else if (layer->type == LAYER_GAUSSIAN)
    {
      tmp1 = (zpos - z) / yanera->parameters.p[layer->sigz];
      sld += tmpsld * exp(-1.0 * tmp1 * tmp1);
    }
    else 
    {
      yanera_error("Unknown component type.", __FILE__, __LINE__, yanera);
    }
    
    layer = layer->next;
  }
  
  /*=========================
   * The last layer 
   * The simple form is a backwards facing error function
   * rho*(1+erf((z-cntr)/(sqrt(2)*sigma)))
   */
  if (layer->isld > -1) tmpsld = yanera->parameters.p[layer->isld];
  else tmpsld = 0.0;
  zpos = yanera->parameters.p[layer->cntr];
  if (yanera->parameters.p[layer->sigz] == 0) 
  {
    if (z >= zpos) sld += tmpsld;
  }
  else
  {
    tmp1 = (z-zpos)/(M_SQRT2*fabs(yanera->parameters.p[layer->sigz]));
    sld  += tmpsld * (1.0 + gsl_sf_erf(tmp1)) / 2.0;
  }    
  return sld;
}
/*============================================================================*/
double functionProfileImag(double z, yanera_model *model, \
                           yanera_container *yanera)
{
  double        thik, sld, tmp1, rho1, rho2;
  yanera_layer *layer1, *layer2;

  /* ------------------------------------
   * The first layer is the layer the neutrons enter from.
   */
  layer1 = model->first_layer;
  if (layer1->isld > -1) sld = yanera->parameters.p[layer1->isld];
  else sld = 0.0;
    
  /* ------------------------------------
   * Loop over the layers.
   */
  layer2 = layer1;
  layer1 = layer1->next;
  thik = 0.0;
  /* ---- */
  while (layer1 != NULL)
  {
    rho1 = 0;
    if (layer1->type == LAYER_FUNCTION)
    {
      rho1 = 0.0;
    }
    else
    {
      if (layer1->isld > -1) rho1 = yanera->parameters.p[layer1->isld];
      else rho1 = 0.0;
      
    }
    
    rho2 = 0;
    if (layer2->type == LAYER_FUNCTION)
    {
      rho2 = 0.0;
    }
    else
    {
      if (layer1->isld > -1) rho2 = yanera->parameters.p[layer1->isld];
      else rho2 = 0.0;
    }    

    if (yanera->parameters.p[layer1->sigz] == 0)
    {
      if (z >= thik) sld += (rho1 - rho2);
    }
    else
    {
      tmp1 = (z-thik)/(M_SQRT2*fabs(yanera->parameters.p[layer1->sigz]));
      sld += (rho1 - rho2) * (1.0 + gsl_sf_erf(tmp1)) / 2.0;
    }

    /* ------------------------------------ 
     * This *if* is because the last layer does not have a thik parameter
     * and without it, we would be accessing the parameter array 
     * outside its bounds.
     */
    if (layer1->thik > -1)
      thik += yanera->parameters.p[layer1->thik];
    
    layer1 = layer1->next;
    layer2 = layer2->next;
  }
  
  return sld;
}
/*============================================================================*/
void minMaxProfile(double *min, double *max, yanera_model *model, \
                   yanera_container *yanera)
{
  double        cntr, sigz;
  yanera_layer *layer;
  
  /* ----------------------------------------
   * Minimum
   */
  layer = model->first_layer;
  min[0] = -1.0;
  if ((yanera->type == MODEL_LAYER) || \
      (yanera->type == MODEL_SLAB) || \
      (yanera->type == MODEL_FUNCTION))
  {
    /*
     *  Use the roughness of the first layer atop the bulk.
     *  MUST BE PRESENT.
     */
    if (yanera->parameters.p[layer->next->sigz] != 0.0)
    { min[0] = -5.0 * fabs(yanera->parameters.p[layer->next->sigz]); }
    else { min[0] -= 1.0; }

  }
  else if (yanera->type == MODEL_COMPONENT)
  { 
    /*
     * Hopefully, the bulk has sigz specified.
     * Start with that. min = cntr - 5*sigz;
     */
    cntr = 0.0;
    if (yanera->parameters.p[layer->sigz] == 0.0) sigz = -1.0;
    else  sigz = -5.0 * fabs(yanera->parameters.p[layer->sigz]);
    min[0] = cntr+sigz;
    /*
     *  Loop over layers looking for objects whose combined
     * center and widths are less than min.
     */
    layer = layer->next;
    while (layer != NULL)
    {
      cntr = yanera->parameters.p[layer->cntr];
      if (layer->type == LAYER_BOX) 
      { cntr -= fabs(yanera->parameters.p[layer->thik]) / 2.0; }
      if (yanera->parameters.p[layer->sigz] == 0.0) sigz = -1.0;
      else  sigz = -5.0 * fabs(yanera->parameters.p[layer->sigz]);
      if (min[0] > (cntr+sigz)) min[0] = cntr + sigz;
      layer = layer->next;
    }
  }
  else 
  {
    yanera_error("Unknown model type in minMaxProfile.",\
                 __FILE__, __LINE__, yanera);
  }
  
  /* ----------------------------------------
   * Maximum
   */
  max[0] = 0.0;
  layer = model->first_layer;
  layer = layer->next;        /* Start one layer down from bulk */
  while (layer->next != NULL) /* Stop on bulk */
  {
    /* Add all thicknesses for slabs, functions */
    if ((yanera->type == MODEL_LAYER) || \
        (yanera->type == MODEL_SLAB)  || \
        (yanera->type == MODEL_FUNCTION))
    {
      max[0] += yanera->parameters.p[layer->thik];
    }
    /*
     *  Look for objects whose combined center and widths are more than max.
     */
    if (yanera->type == MODEL_COMPONENT)
    {
      cntr = yanera->parameters.p[layer->cntr];
      if (layer->type == LAYER_BOX) 
      { cntr += fabs(yanera->parameters.p[layer->thik]) / 2.0; }
      if (yanera->parameters.p[layer->sigz] == 0.0) sigz = 1.0;
      else  sigz = 5.0 * fabs(yanera->parameters.p[layer->sigz]);
      if (max[0] < (cntr+sigz)) max[0] = cntr + sigz;
    }
    layer = layer->next;
  }
  /* We're now at the last layer, a bulk layer */
  /* Add width of transition to last layer */
  if ((yanera->type == MODEL_LAYER) || \
      (yanera->type == MODEL_SLAB)  || \
      (yanera->type == MODEL_FUNCTION))
  {
    if (yanera->parameters.p[layer->sigz] == 0.0) max[0] += 1.0;
    else max[0] += 5.0 * yanera->parameters.p[layer->sigz];
  }
  /* Is the combined center and width of the last layer more than max? */
  if (yanera->type == MODEL_COMPONENT)
  {
    cntr = yanera->parameters.p[layer->cntr];
    if (layer->type == LAYER_BOX) 
    { cntr += fabs(yanera->parameters.p[layer->thik]) / 2.0; }
    if (yanera->parameters.p[layer->sigz] == 0.0) sigz = 1.0;
    else  sigz = 5.0 * fabs(yanera->parameters.p[layer->sigz]);
    if (max[0] < (cntr+sigz)) max[0] = cntr + sigz;
  }

  return;
}

---