/*ufunct one of K. Volk's quantum heating routines for grains */
#include "cddefines.h"
#include "grainvar.h"
#include "lastdt.h"
#include "ufunct.h"

double ufunct(double temp, 
  long int ndust1)
{
	long int itype1, 
	  j;
	double ufunct_v;
	static double power[4]={3.0,2.3,1.68,1.0};
	static double tlim[5]={0.,50.,150.,500.,1.e30};
	static double cval[4]={2.300e-10,4.7141e-09,1.4081e-07,1.6806e-05};
	static double term[4]={0.,2.8749e-05,4.6752e-04,7.4916e-03};
	static double atoms[NDUST]={3.991e07,1.725e07,5.863e08,2.534e08,
	  3.902e05,3.902e07,1.725e07,1.725e07,24.0,144.0,0.0,0.0,0.0,0.0,
	  0.0,0.0,0.0,0.0,0.0,0.0};

#	ifdef DEBUG_FUN
	fputs( "<+>ufunct()\n", debug_fp );
#	endif

	/* this common block holds the number of atoms and base grain type
	 * for the 10 grains that can be read in from files. */
	/* coefficients converted to rydberg energy units, per unit atom
	 * assuming a density of 3.3 g/cm^3 and pure silicon dioxide.
	 * this is not right, but the result is correct because atoms 
	 * was calculated using the same assumption. */

	if( ndust1 > NDUST - 10 )
	{
		atoms[ndust1-1] = lastdt.anumkv[ndust1-11];
		itype1 = lastdt.itype[ndust1-11] + 1;
	}
	else
	{
		itype1 = 0;
	}

	if( ((ndust1 == 1 || ndust1 == 3) || ndust1 > 7) || itype1 == 
	  1 )
	{
		/* pah/graphitic grains (and the greybody grains) */
		ufunct_v = (pow(temp,3.3))*1.904e-11/(1. + 6.51e-03*temp + 
		  1.5e-06*temp*temp + 8.3e-07*(pow(temp,2.3)));
	}
	else
	{
		/* silicate grains
		 * do j=1,4 */
		j = 1;
		ufunct_v = 0.;
		while( j < 4 && ufunct_v == 0 )
		{
			/* limits of tlim set above, 0 and 1e30, so must happen */
			if( temp > tlim[j-1] && temp <= tlim[j] )
			{
				ufunct_v = term[j-1] + cval[j-1]*(pow(temp,power[j-1]) - 
				  pow(tlim[j-1],power[j-1]));
				/* go to 455 */
			}
			j += 1;
		}
	}
	/*455  ufunct=ufunct*atoms(ndust1) */
	ufunct_v *= atoms[ndust1-1];


#	ifdef DEBUG_FUN
	fputs( " <->ufunct()\n", debug_fp );
#	endif
	return( ufunct_v );
}