/*ParseGrain parse parameters on grains command */
#include <string.h>
#include <limits.h>
#include "cddefines.h"
#include "grainvar.h"
#include "physok.h"
#include "input.h"
#include "dheton.h"
#include "ffmtread.h"
#include "lgmatch.h"
#include "getquote.h"
#include "parse.h"
#include "lastdt.h"
/*rdfile read grains parameters from a file. called by getgrains when grain number  >10 */
static void rdfile( long int igr , /* integer counter for grain type we will read in */
	char chFile[] , /* name of file with data */
	int lgBinary/* is this a binary file?  ascii if false */
	); 

void ParseGrain(char *chCard, 
  int *lgDset)
{
	int lgEOL, 
	  lgEOL1, 
	  lgLinSet, 
	  lgLogOn, 
	  lgVaryAbun;

	char chFile[200];/*possible name of input file with opacities */
	long int i, 
	  ipnt;
	double dep, 
	  dep1;

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

	/* grain model, numbers are logs of grain abundance rel to std */
	GrainVar.lgDustOn = TRUE;
	*lgDset = TRUE;

	/* log - linear option for grains */
	lgLogOn = FALSE;
	lgLinSet = FALSE;
	if( lgMatch(" LOG",chCard) )
	{
		lgLogOn = TRUE;
		lgLinSet = FALSE;
	}

	else if( lgMatch("LINE",chCard) )
	{
		lgLogOn = TRUE;
		lgLinSet = TRUE;
	}

	/* check for the keyword "function" - this is a special case */
	if( lgMatch("FUNC",chCard) )
	{
		lgVaryAbun = TRUE;
	}
	else
	{
		lgVaryAbun = FALSE;
	}

	/* get the grain abundance as the first parameter */
	i = 5;
	dep = FFmtRead(chCard,&i,76,&lgEOL);

	/* was keyword log on the line? */
	if( lgLogOn )
	{
		/* log or linear was specified, which was it */
		if( lgLinSet )
		{
			/* linear quantity enterd, and we want a log */
			if( dep > 0. )
			{
				dep = log10(dep);
			}
			else
			{
				fprintf( ioQQQ, " Impossible for linear negative abundance.\n" );
				fprintf( ioQQQ, " Abundance entered was%10.2e\n", dep );
				fprintf( ioQQQ, " Sorry.\n" );
				puts( "[Stop in getgrain]" );
				exit(1);
			}
		}
	}

	else
	{
		/* neither log nor linear specified, check sign
		 * force it to be a log - linear if greater than 0 */
		if( dep > 0. )
		{
			dep = log10(dep);
		}
	}

	/* now change to linear since will be used throughout code */
	dep = pow(10.,dep);

	/* it is possible that there is a file name on the command line - 
	 * if so then we want to call correct reoutine, and not look for keywords 
	 * the signature of a eyword is a pair of quotes - is one present?
	 */
	if( strchr( input.chOrgCard,'\"' ) != NULL )
	{
		/* 
		 * this will both scan in whatever label is inside the quotes in OrgCard, 
		 * but also remove the contents there and in chCard,
		 * so that following keywords will not trigger off it
		 */
		GetQuote( chFile , chCard );
	}
	else
	{
		/* set end of line sentinal so we know file not on line */
		chFile[0] = '\0';
	}

	/* 
	 * NB! lgDustOn1(nd) turns on grains for some species.  None of these
	 * commands turn off other grains, so successive commands will turn on
	 * more and more grains.  This is intended
	 * 
	 */
	if( lgMatch("ORIO",chCard) )
	{
		/* optional keyword ORION to use orion curves for large R grains
		 * >>chng 96 nov 30, into loop for extra K Volk grain types */
		GrainVar.lgDustOn1[2] = TRUE;
		GrainVar.lgDustOn1[3] = TRUE;
		GrainVar.dstfactor[2] = (float)dep;
		GrainVar.dstfactor[3] = (float)dep;
		/* option to have grain abundance depend on depth */
		if( lgVaryAbun )
		{
			GrainVar.lgDustVary[2] = TRUE;
			GrainVar.lgDustVary[3] = TRUE;
		}
	}

	else if( lgMatch(" PAH",chCard) )
	{
		/* >>chng 97 mar 1, added this keyworkd
		 * optional keyword  PAH to use two PAH grain species */
		GrainVar.lgDustOn1[8] = TRUE;
		GrainVar.lgDustOn1[9] = TRUE;
		GrainVar.dstfactor[8] = (float)dep;
		GrainVar.dstfactor[9] = (float)dep;

		/* for PAH grains (only) there is a qheat option for quantum heating */
		if( lgMatch("QHEA",chCard) )
		{
			GrainVar.lgQHeat[8] = TRUE;
			GrainVar.lgQHeat[9] = TRUE;
		}

		/* option to have grain abundance depend on depth */
		if( lgVaryAbun )
		{
			GrainVar.lgDustVary[8] = TRUE;
			GrainVar.lgDustVary[9] = TRUE;
		}
	}

	else if( lgMatch("GRAY",chCard) || lgMatch("GREY",chCard) )
	{
		/* >>chng 97 mar 1, added this keyworkd, 97 july 18, second spelling
		 * optional keyword GRAY to use gray grain, 0.1 micron */
		GrainVar.lgDustOn1[7] = TRUE;
		GrainVar.dstfactor[7] = (float)dep;

		/* option to have grain abundance depend on depth */
		if( lgVaryAbun )
		{
			GrainVar.lgDustVary[7] = TRUE;
		}
	}

	else if( lgMatch(" AGB",chCard) || lgMatch("AGB ",chCard) )
	{
		/* key used to be planetary nebula, should be agb
		 * planetary nebula - post agb, from Kevin Volk
		 * Try to read the second depletion factor (K Volk) */
		dep1 = FFmtRead(chCard,&i,76,&lgEOL1);
		if( lgEOL1 )
		{
			dep1 = dep;
		}

		/* >>chng 97 july 17, linear dep1
		 * if(dep1.gt.0.) dep1 = LOG10(dep1) */
		if( dep1 <= 0. )
		{
			dep1 = pow(10.,dep1);
		}

		/* End of changes here (K Volk) */
		GrainVar.lgDustOn1[0] = TRUE;
		GrainVar.lgDustOn1[6] = TRUE;
		GrainVar.dstfactor[0] = (float)dep;

		/* >>chng 96 nov 30, as per K Volk second parameter */
		GrainVar.dstfactor[6] = (float)dep1;

		/* option to have grain abundance depend on depth */
		if( lgVaryAbun )
		{
			GrainVar.lgDustVary[0] = TRUE;
			GrainVar.lgDustVary[6] = TRUE;
		}
	}

	else
	{
		/* turn one type of dust on if second number present */
		ipnt = (long int)FFmtRead(chCard,&i,76,&lgEOL);
		/* leave documented behavior or grain pointers as they were,
		 * so this is one greater than array element.  decrement here 
		 * to bring onto C scale */
		--ipnt;

		/* lgEOL if no second number was present */
		if( lgEOL )
		{
			/* no second number, use default ISM dust with optional depletion
			 * Mathis et al. ISM grains */
			GrainVar.lgDustOn1[0] = TRUE;
			GrainVar.lgDustOn1[1] = TRUE;
			GrainVar.dstfactor[0] = (float)dep;
			GrainVar.dstfactor[1] = (float)dep;

			/* option to have grain abundance depend on depth */
			if( lgVaryAbun )
			{
				GrainVar.lgDustVary[0] = TRUE;
				GrainVar.lgDustVary[1] = TRUE;
			}
		}

		/* NDUST is 20 in grainvar.h */
		else if( ipnt < 0 || ipnt >= NDUST )
		{
			fprintf( ioQQQ, " Pointer to grain type bad.  Sorry.\n" );
			puts( "[Stop in getgrain]" );
			exit(1);
		}

		else
		{
			/* this branch if grain number on line */
			GrainVar.lgDustOn1[ipnt] = TRUE;
			GrainVar.dstfactor[ipnt] = (float)dep;

			/* option to have grain abundance depend on depth */
			if( lgVaryAbun )
			{
				GrainVar.lgDustVary[ipnt] = TRUE;
			}

			/* Changes to dust for other grains (K Volk)
			 * this branch if grain species greater than 10 - we will
			 * read in optical properties from ancillary file */
			if( ipnt >= 11 )
			{
				int lgBinary;
				if( lgMatch("BINA",chCard) )
				{
					lgBinary = TRUE;
				}
				else
				{
					lgBinary = FALSE;
				}
				rdfile( ipnt , chFile , lgBinary );
			}

			/* reads in parameters for the grain from a file if the
			 * number is >= 11 */

			/* Add change for quantum heating, keeping track of depletions of grain */
			GrainVar.qdep[ipnt] *= (float)pow(10.,dep);

			/* Add qheat option
			 * for grains set by number (only) there is a qheat option for
			 * quantum heating */
			if( lgMatch("QHEA",chCard) )
			{
				GrainVar.lgQHeat[ipnt] = TRUE;
			}

			/* End of changes here (K Volk) */
		}
	}

	/* option to turn off photoelectric heating by grain, NO HEATING */
	if( lgMatch("O HE",chCard) )
	{
		physok.lgPhysOK = FALSE;
		dheton.lgDHetOn = FALSE;
	}

	/* option to turn off gas cooling by grain, NO COOLING */
	if( lgMatch("O CO",chCard) )
	{
		physok.lgPhysOK = FALSE;
		dheton.lgDColOn = FALSE;
	}

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

/*rdfile read grains parameters from a file. called by getgrain when grain number  >10 */

static void rdfile( long int igr,
				 char chLine[] ,/* the name of the file */
				 int lgBinary /* is it a binary file?*/
				 )
{
	FILE * ioFILE ;
	long int  i,
		nfield ;
	float s1;
	char chCard[133];

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

	/* read grains parameters from a file.  It is called when the grain number  >11
	 * The sublimation temperature is a new addition for c90.
	 * This subroutine added by (K Volk) 27 August 1993 */

	/* first try to open the file, 
	 * its name was read in with GetQuote in GetGrain, but only it quotes were entered */

	/* if first char is null then we know name was not properly entered */
	if( chLine[0]=='\0' )
	{
		/* big problem - we got here but double quotes were not picked up -
		 * not enough information to proceed */
		fprintf( ioQQQ , " rdfile got grain species that needed input file, but none specified.\n");
		fprintf( ioQQQ , " file name must occur between pair of double quotes.\n");
		puts( "[Stop in rdfile]" );
		exit(1);
	}
	if( lgBinary )
	{
typedef struct {
	double re;
	double im;
} complex;

/* maximum no. of parameters for grain size distribution */
#define NSD 5

typedef struct {
	double a[NSD];  /* parameters for size distribution */
	double *rr;     /* rr[nn]: abcissas for Gauss quadrature */
	double *ww;     /* ww[nn]: weights for Gauss quadrature */
	double area;    /* average grain surface area (darea) */
	double vol;     /* average grain volume (dsize = 3.*vol/area) */
	long int nn;    /* no. of abcissas used in Gauss quadrature (per partition) */
	long int nPart; /* no. of partitions for size distribution */
	int sdCase;     /* 0 = single size, 1 = simple powerlaw size distr, ... */
	int lgLogScale; /* use logarithmic mesh for integration over size ? */
} sd_data;

/* maximum no. of principal axes for crystalline grains */
#define NAX 3

typedef struct {
	double *wavlen[NAX];
	complex *n[NAX];
	double *nr1[NAX];
	double wt[NAX];
	double abun;            /* dustp(3)*dustp(4) */
	double depl;            /* dustp(4) */
	double elmAbun[LIMELM]; /* dstela */
	double mol_weight;      /* dustp(2) */
	double rho;             /* dustp(1) */
	double work;            /* DustWorkFcn */
	double subl_temp;       /* sublimat */
	long int ndata[NAX];
	long int nAxes;
} grain_data;

	int lgErr;
	long int n,
	  nfreqp;
	double anu[NCELL],
	  abs_cs[NCELL],
	  sct_cs[NCELL];
	size_t num;
	char *fnam;
	sd_data sd;
	grain_data gd;
	FILE *ioGRAIN;

	/* open the data file containing the compiled opacity data */
	fnam = "d:\\projects\\cloudy\\mie\\test.po";
	if( (ioGRAIN = fopen(fnam, "rb")) == NULL ) 
	{
		fprintf( ioQQQ, "error opening file: %s\n",fnam );
		puts( "[Stop in read grains]" );
		exit(1);
	}

	lgErr = FALSE;
	num = 1;
	/* sd is size distribution data structure */
	lgErr = lgErr || ( fread(&sd,sizeof(sd_data),num,ioGRAIN) != num );

	/* gd is grain data structure */
	lgErr = lgErr || ( fread(&gd,sizeof(grain_data),num,ioGRAIN) != num );

	/* n is number of energy array as written out, was == NCELL in version of
	 * code that wrote it, should still be NCELL here */
	n = INT_MIN;
	lgErr = lgErr || ( fread(&n,sizeof(long),num,ioGRAIN) != num );
	if( !lgErr && n != NCELL ) 
	{
		fprintf( ioQQQ, " Incompatible array length in file: %s\n",fnam );
		fprintf( ioQQQ, " Please recompile this file\n" );
		fclose(ioGRAIN);
		puts( "[Stop]" );
		exit(1);
	}

	/* nfreqp is no. of cells that are actually initialized */
	lgErr = lgErr || ( fread(&nfreqp,sizeof(long),num,ioGRAIN) != num ); 
	num = NCELL;
	lgErr = lgErr || ( fread(anu,sizeof(double),num,ioGRAIN) != num );
	lgErr = lgErr || ( fread(abs_cs,sizeof(double),num,ioGRAIN) != num );
	lgErr = lgErr || ( fread(sct_cs,sizeof(double),num,ioGRAIN) != num );
	fclose(ioGRAIN);
	if( lgErr ) 
	{
		fprintf( ioQQQ, " Error reading file: %s\n", fnam );
		fprintf( ioQQQ, " Please recompile this file\n" );
		puts( "[Stop]" );
		exit(1);
	}
	strcpy( chGrainVar.chDstLab[igr] , "mie " );
	GrainVar.DustWorkFcn[igr] = (float)gd.work;
	/* 0 for graphite, 1 for silicate */
	lastdt.itype[igr-10] = 0;
	/* number of atoms in the grain */
	lastdt.anumkv[igr-10] = 0.;
	GrainVar.qhnrat[igr] = 0.;
	GrainVar.Tsublimat[igr] = (float)gd.subl_temp;

	GrainVar.dustp[igr][0] = (float)gd.rho;
	GrainVar.dustp[igr][1] = (float)gd.mol_weight;
	GrainVar.dustp[igr][3] = (float)gd.depl;
	GrainVar.dustp[igr][2] = (float)(gd.abun/gd.depl);
	/* this shows that exact values are read in for the energy array in place */
	GrainVar.ndpts[igr] = -1;
	GrainVar.darea[igr] = (float)sd.area;
	GrainVar.dsize[igr] = (float)(3.*sd.vol/sd.area);

	/* now copy the scattering and absorption arrays over */
	for( i=0; i<nfreqp; ++i )
	{
		GrainVar.dstab1[igr][i] = abs_cs[i];
		GrainVar.dstsc1[igr][i] = sct_cs[i];
		{
			enum {DEBUG=FALSE};
			if( DEBUG )
			{
				fprintf( ioQQQ,"%.2e %.2e %.2e\n",
					anu[i], GrainVar.dstab1[igr][i] ,GrainVar.dstsc1[igr][i] );
			}
		}
	}

	/* now copy the scattering and absorption arrays over */
	for( i=nfreqp; i<NCELL; ++i )
	{
		GrainVar.dstab1[igr][i] = 0.;
		GrainVar.dstsc1[igr][i] = 0.;
	}
	}
	else
	{

	ioFILE = fopen( chLine , "r" );
	if( ioFILE == NULL )
	{
		fprintf( ioQQQ , "rdfile could not open grain opacity file named %s\n",chLine);
		puts( "[Stop in rdfile]" );
		exit(1);
	}

	s1 = GrainVar.Tsublimat[igr];

	/* read first line */
	if( fgets( chCard , sizeof(chCard) , ioFILE ) == NULL )
	{
		fprintf( ioQQQ, " rdfile error getting first line\n" );
		puts( "[Stop in rdfile]" );
		exit(1);
	}

	nfield = sscanf( chCard ,  " %10s%f %1li %12e %12e %7f",
	  chGrainVar.chDstLab[igr], 
	  &GrainVar.DustWorkFcn[igr], 
	  &lastdt.itype[igr-10], 
	  &lastdt.anumkv[igr-10], 
	  &GrainVar.qhnrat[igr], 
	  &GrainVar.Tsublimat[igr] );

	/* mak sure we got 6 values  */
	if( nfield != 6 )
	{
		fprintf( ioQQQ, " rdfile error parsing first line\n" );
		puts( "[Stop in rdfile]" );
		exit(1);
	}

	/* The "old" dust input file does not have the sublimation temperature 
	 * so it is reset to the initial value if it is zero, or if the value is 
	 * absurdly small.
	 */
	if( GrainVar.Tsublimat[igr] <= 50. )
	{
		GrainVar.Tsublimat[igr] = s1;
	}

	/* read second line */
	if( fgets( chCard , sizeof(chCard) , ioFILE ) == NULL )
	{
		fprintf( ioQQQ, " rdfile error getting second line\n" );
		puts( "[Stop in rdfile]" );
		exit(1);
	}

	nfield = sscanf( chCard ,"%9e %9e %9e  %9e %3li  %9e  %9e",
		&GrainVar.dustp[igr][0] ,
		&GrainVar.dustp[igr][1] ,
		&GrainVar.dustp[igr][2] ,
		&GrainVar.dustp[igr][3] ,
		&GrainVar.ndpts[igr], 
		&GrainVar.darea[igr], 
	   &GrainVar.dsize[igr]
		);
	/* make sure we got 7 values */
	if( nfield != 7 )
	{
		fprintf( ioQQQ, " rdfile error parsing second line\n" );
		puts( "[Stop in rdfile]" );
		exit(1);
	}

	if( GrainVar.ndpts[igr] > LIMCRS )
		GrainVar.ndpts[igr] = LIMCRS;

	/* read third line */
	if( fgets( chCard , sizeof(chCard) , ioFILE ) == NULL )
	{
		fprintf( ioQQQ, " rdfile error getting third line\n" );
		puts( "[Stop in rdfile]" );
		exit(1);
	}

/*		original fortran
      read(89,'(4(f5.1,1x,e12.4))', err=585,end=585) (dstelz(j,igr),
     1                              dstela(j,igr),j=1,4)
	  */
	nfield = sscanf( chCard ,"%5f %12e%5f %12e%5f %12e%5f %12e",
		&GrainVar.dstelz[igr][0], 
		&GrainVar.dstela[igr][0],
		&GrainVar.dstelz[igr][1], 
		&GrainVar.dstela[igr][1],
		&GrainVar.dstelz[igr][2], 
		&GrainVar.dstela[igr][2],
		&GrainVar.dstelz[igr][3], 
		&GrainVar.dstela[igr][3]
		);
	/* make sure we got 8 values */
	if( nfield != 8 )
	{
		fprintf( ioQQQ, " rdfile error parsing third line\n" );
		puts( "[Stop in rdfile]" );
		exit(1);
	}

	for( i=0; i < GrainVar.ndpts[igr]; i++ )
	{
		/* read the data */
		if( fgets( chCard , sizeof(chCard) , ioFILE ) == NULL )
		{
			fprintf( ioQQQ, " rdfile error getting data, line=%li\n",i );
			puts( "[Stop in rdfile]" );
			exit(1);
		}

		nfield = sscanf( chCard ,"%11e%11e%11e",
			&GrainVar.eev[igr][i], 
			&GrainVar.sab[igr][i], 
			&GrainVar.sse[igr][i]
			);
		/* make sure we got 8 values */
		if( nfield != 3 )
		{
			fprintf( ioQQQ, " rdfile error parsing data line, %li=\n",i );
			puts( "[Stop in rdfile]" );
			exit(1);
		}
	}

	fclose( ioFILE );
	}

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