/*
 * SUMMARY:      SeparateRadiation.c - Separate radiation
 * USAGE:        Part of DHSVM
 *
 * AUTHOR:       Bart Nijssen
 * ORG:          University of Washington, Department of Civil Engineering
 * E-MAIL:       nijssen@u.washington.edu
 * ORIG-DATE:    Apr-96
 * LAST-MOD: Mon Apr 14 12:59:20 1997 by Bart Nijssen <nijssen@meter.ce.washington.edu>
 * DESCRIPTION:  Separate observed downward solar radiation into diffuse and 
 *               direct beam radiation based on the clearness index (kt)
 * DESCRIP-END.
 * FUNCTIONS:    SeparateRadiation()
 * COMMENTS:     
 */

#ifndef lint
static char vcid[] = "$Id: SeparateRadiation.c,v 1.2 2001/03/31 00:10:35 waichler Exp $";
#endif /* lint */

#include <math.h>
#include "settings.h"
#include "rad.h"

/*****************************************************************************
  Function name: SeparateRadiation()

  Purpose      : Separate observed downward solar radiation into diffuse and 
                 direct beam radiation based on the clearness index (kt)

  Required     :
    float TotalSolar - Total amount of solar (shortwave) radiation (W/m2)
    float ClearIndex - Clearness index
    float *Beam      - Direct beam radiation (W/m2)
    float *Diffuse   - Diffuse radiation (W/m2)

  Returns      : void

  Modifies     :
    float *Beam      - Direct beam radiation (W/m2)
    float *Diffuse   - Diffuse radiation (W/m2)

  Comments     : Based on: Erbs, D.G., S.A. Klein, and J.A. Duffie, 
                 "Estimation of the diffuse fraction for hourly, daily and 
                 monthly-average global radiation", Solar Energy, V.28, n.4, 
                 pp.293-302, 1982
*****************************************************************************/
void SeparateRadiation(float TotalSolar, float ClearIndex, 
		       float *Beam, float *Diffuse)
{   
  /* The following relationships were developed for hourly data, and are here
     used varying timesteps (Time.Dt).  It would probably be better to
     develop explicit relationships for different timesteps */
  if (ClearIndex <= 0.22) {
    *Diffuse = (1 - 0.09 * ClearIndex);
  } 
  else if (ClearIndex <= 0.80) {
    *Diffuse = 0.9511 - 0.1604 * ClearIndex + 
      4.388 * ClearIndex*ClearIndex -
      16.638 * ClearIndex*ClearIndex*ClearIndex + 
      12.366 * ClearIndex*ClearIndex*ClearIndex*ClearIndex;
  } else
    *Diffuse = 0.165;
  
  *Diffuse *= TotalSolar;
  *Beam     = TotalSolar - *Diffuse;
}
