/*PunOpac punch total opacity in any element, punch opacity command */
#include <string.h>
#include "cddefines.h"
#include "physconst.h"
#include "pnunit.h"
#include "nhe1lvl.h"
#include "opacpoint.h"
#include "optype.h"
#include "rfield.h"
#include "kshllenr.h"
#include "opac.h"
#include "phycon.h"
#include "he.h"
#include "nhe1.h"
#include "nh.h"
#include "hydrogenic.h"
#include "linelabl.h"
#include "heavy.h"
#include "ionfracs.h"
#include "elementnames.h"
#include "opacz.h"
#include "punopac.h"
#include "printe82.h"
#include "opinb.h"
#include "putopacity.h"

void PunOpac(FILE * io, 
  long int np)
{
	/* this will be the file name for opacity output */
	char chFileName[100];	

	/* this is its pointer */
	FILE *ioFileName;			

	char chNumbers[31][3] = {"0","1","2","3","4","5","6","7","8","9",
		"10","11","12","13","14","15","16","17","18","19",
		"20","21","22","23","24","25","26","27","28","29","30"};

	long int i, 
	  ilow, 
	  ion, 
	  ipS, 
	  ipZ, 
	  j, 
	  nelem;

	double ener, 
	  ener3;

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

	/* punch total opacity in any element, punch opacity command
	 * io is io unit number, np is pointer within stack of punches */

	if( strcmp(optype.chOpcTyp,"TOTL") == 0 )
	/* total opacity */
	{
		for( j=1; j <= rfield.nflux; j++ )
		{
			fprintf( io, "%12.4e%10.2e%10.2e%10.2e%10.2e %4.4s\n", 
			  rfield.anu[j-1], 
			  opac.opac[j-1] + opac.scatop[j-1], 
			  opac.opac[j-1], 
			  opac.scatop[j-1], 
			  opac.scatop[j-1]/
			  MAX2(1e-37,opac.scatop[j-1]+opac.opac[j-1]), 
			  LineLabl.chContLabel[j-1] );
		}

		fprintf( io, "\n" );
	}

	/* subshell photo cross sections */
	else if( strcmp(optype.chOpcTyp,"SHEL") == 0 )
	{
		nelem = (long)pnunit.punarg[0][np];
		ion = (long)pnunit.punarg[1][np];
		ipS = (long)pnunit.punarg[2][np];
		for( i=OpacPoint.ipElement[0][ipS-1][ion-1][nelem-1]-1; 
		  i < OpacPoint.ipElement[1][ipS-1][ion-1][nelem-1]; i++ )
		{
			fprintf( io, 
				"%11.3e %11.3e\n", rfield.anu[i], 
			  opac.OpacStack[i-OpacPoint.ipElement[0][ipS-1][ion-1][nelem-1]+
			  OpacPoint.ipElement[2][ipS-1][ion-1][nelem-1]] );
		}
	}

	/* figure for hazy */
	else if( strcmp(optype.chOpcTyp,"FIGU") == 0 )
	{
		ipZ = 0;
		for( i=nh.ipHn[0][0]-1; i < (nhe1Com.nhe1[0] - 1); i++ )
		{
			ener = rfield.anu[i]*0.01356;
			ener3 = 1e24*POW3(ener);
			fprintf( io, 
				"%12.4e%12.4e%12.4e%12.4e%12.4e\n", 
			  rfield.anu[i], ener, 
			  opac.OpacStack[i-nh.ipHn[0][IP1S]+ OpacPoint.ipHOpac[ipZ][IP1S]]*ener3, 
			  0., 
			  opac.opac[i]/ phycon.hden*ener3 );
		}

		for( i=nhe1Com.nhe1[0]-1; i < rfield.nupper; i++ )
		{
			ener = rfield.anu[i]*0.01356;
			ener3 = 1e24*POW3(ener);
			fprintf( io, 
				"%12.4e%12.4e%12.4e%12.4e%12.4e\n", 
			  rfield.anu[i], 
			  ener, 
			  opac.OpacStack[i-nh.ipHn[0][IP1S]+ OpacPoint.ipHOpac[ipZ][IP1S]]*ener3, 
			  opac.OpacStack[i-nhe1Com.nhe1[0]+ OpacPoint.iophe1[0]]*he.ahe/phycon.hden*ener3, 
			  opac.opac[i]/phycon.hden*ener3 );
		}
	}

	/* hydrogen */
	else if( strcmp(optype.chOpcTyp,"HYDR") == 0 )
	{
		ipZ = 0;
		opacz();
		opinb(OpacPoint.ipHOpac[ipZ][IP1S],nh.ipHn[ipZ][IP1S],
		  rfield.nupper,1.,0.,'v');
		ilow = Heavy.ipHeavy[0][0];

		/* start filename for output */
		strcpy( chFileName , elementnames.chElementNameShort[0] );

		/* continue filename with ionization stage */
		strcat( chFileName , chNumbers[1] );

		/* end it wil .opc, name will be of form HYDR.opc */
		strcat( chFileName , ".opc" );

		/* now try to open it for writing */
		if( NULL==(ioFileName = fopen( chFileName , "w" )) )
		{
			fprintf( ioQQQ," could not open opacity file %s for writing.\n",chFileName);
			puts( "[Stop in punopac]" );
			exit(1);
		}

		for( j=ilow; j <= rfield.nupper; j++ )
		{
			/* photon energy in rydbergs */
			PrintE93( ioFileName , rfield.anu[j-1]*EVRYD );
			fprintf( ioFileName , " ");
			/* cross section in megabarns */
			PrintE93( ioFileName, opac.opac[j-1]/1e-18 );
			fprintf( ioFileName , "\n");
		}

		fclose( ioFileName );
		puts( "[Stop in punopac]" );
		exit(1);

	}

	/* helium */
	else if( strcmp(optype.chOpcTyp,"HELI") == 0 )
	{
		/* atomic helium first, HELI1.opc */
		opacz();
		opinb(OpacPoint.iophe1[0],nhe1Com.nhe1[0],rfield.nflux,1.,0.,'v');
		ilow = Heavy.ipHeavy[0][1];

		/* start filename for output */
		strcpy( chFileName , elementnames.chElementNameShort[1] );

		/* continue filename with ionization stage */
		strcat( chFileName , chNumbers[1] );

		/* end it wil .opc, name will be of form HYDR.opc */
		strcat( chFileName , ".opc" );

		/* now try to open it for writing */
		if( NULL==(ioFileName = fopen( chFileName , "w" )) )
		{
			fprintf( ioQQQ," could not open opacity file %s for writing.\n",chFileName);
			puts( "[Stop in punopac]" );
			exit(1);
		}

		for( j=ilow; j <= rfield.nupper; j++ )
		{
			/* photon energy in rydbergs */
			PrintE93( ioFileName , rfield.anu[j-1]*EVRYD );
			fprintf( ioFileName , " ");
			/* cross section in megabarns */
			PrintE93( ioFileName, opac.opac[j-1]/1e-18 );
			fprintf( ioFileName , "\n");
		}
		fclose( ioFileName );

		/* now do helium ion, HELI2.opc */
		opacz();
		opinb(OpacPoint.ipHOpac[1][IP1S],nh.ipHn[1][IP1S],rfield.nupper,1.,0.,'v');
		ilow = Heavy.ipHeavy[1][1];

		/* start filename for output */
		strcpy( chFileName , elementnames.chElementNameShort[1] );

		/* continue filename with ionization stage */
		strcat( chFileName , chNumbers[2] );

		/* end it wil .opc, name will be of form HYDR.opc */
		strcat( chFileName , ".opc" );

		/* now try to open it for writing */
		if( NULL==(ioFileName = fopen( chFileName , "w" )) )
		{
			fprintf( ioQQQ," could not open opacity file %s for writing.\n",chFileName);
			puts( "[Stop in punopac]" );
			exit(1);
		}

		for( j=ilow; j <= rfield.nupper; j++ )
		{
			/* photon energy in rydbergs */
			PrintE93( ioFileName , rfield.anu[j-1]*EVRYD );
			fprintf( ioFileName , " ");
			/* cross section in megabarns */
			PrintE93( ioFileName, opac.opac[j-1]/1e-18 );
			fprintf( ioFileName , "\n");
		}

		fclose( ioFileName );
		puts( "[Stop in punopac]" );
		exit(1);
	}

	else
	{
		/* check for hydroge through zinc, nelem is atomic number on the c scale */
		nelem = -1;
		i = 0;
		while( i < LIMELM )
		{
			if( strcmp(optype.chOpcTyp,elementnames.chElementNameShort[i]) == 0 )
			{
				nelem = i;
				break;
			}
			++i;
		}

		/* nelem is still negative if above loop fell through */
		if( nelem < 0 )
		{
			fprintf( ioQQQ, " Unidentified opacity key=%4.4s\n", 
			  optype.chOpcTyp );
			puts( "[Stop in punopac]" );
			exit(1);
		}

		/* generic driving of PutOpacity */
		hydro.hn[nelem][IP1S] = 1.;
		hydro.hbn[nelem][IP1S] = 0.;

		for( i=0; i <= nelem; i++ )
		{
			for( j=0; j < (nelem + 2); j++ )
			{
				xIonFracs[j][nelem] = 0.;
			}

			xIonFracs[i+1][nelem] = 1.;

			opacz();

			/* generate opacity with standard routine - this is the one
			 * called in addopac to make opacities in usual calculations */
			PutOpacity(nelem);

			/* start filename for output */
			strcpy( chFileName , elementnames.chElementNameShort[nelem] );

			/* continue filename with ionization stage */
			strcat( chFileName , chNumbers[i+1] );

			/* end it wil .opc, name will be of form HYDR.opc */
			strcat( chFileName , ".opc" );

			/* now try to open it for writing */
			if( NULL==(ioFileName = fopen( chFileName , "w" )) )
			{
				fprintf( ioQQQ," could not open opacity file %s for writing.\n",chFileName);
				puts( "[Stop in punopac]" );
				exit(1);
			}

			ilow = Heavy.ipHeavy[i][nelem];

			for( j=ilow-1; j < KshllEnr.KshellLimit; j++ )
			{
				/* photon energy in rydbergs */
				PrintE93( ioFileName , rfield.anu[j]*EVRYD );
				fprintf( ioFileName , " ");

				/* cross section in megabarns */
				PrintE93( ioFileName, opac.opac[j]/1e-18 );
				fprintf( ioFileName , "\n");
			}
			/* close this ionization stage */
			fclose( ioFileName );
		}

		puts( "[Stop in punopac]" );
		exit(1);
	}


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