/* This file is part of Cloudy and is copyright (C)1978-2013 by Gary J. Ferland and
 * others.  For conditions of distribution and use see copyright notice in license.txt */
/*plot master routine to generate some sort of plot */
#include "cddefines.h"
#include "iterations.h"
#include "called.h"
#define	IHI	59
#define	IWID	121
#include "input.h"
#include "rfield.h"
#include "trace.h"
#include "radius.h"
#include "geometry.h"
#include "opacity.h"
#include "dense.h"
#include "hcmap.h"
#include "plot.h"

t_plotCom plotCom;

/*pltcon generate plot of continuum array */
STATIC void pltcon(long int np, 
  const char *chCall);

/*pltmap generate plot of heating and cooling map */
STATIC void pltmap(long int np, 
  const char *chCall);

/*pltopc generate plot of local gas opacity */
STATIC void pltopc(long int np, 
  const char *chCall);

/* this is the base routine that actually makes the plots, called by above */
STATIC void pltr(realnum[],realnum[],long,double,double,double,double,
	char,char*,long,bool);


void plot(const char *chCall)
{
	long int np;

	DEBUG_ENTRY( "plot()" );

	/* return if this is not the last iteration, or a plot not required,
	 * or we are not speaking */
	if( !plotCom.lgPlotON || !called.lgTalk )
	{ 
		return;
	}

	if( !iterations.lgLastIt && (strcmp(chCall,"FIRST") != 0) )
	{ 
		return;
	}

	/* loop over all the requested plots */
	for( np=0; np < plotCom.nplot; np++ )
	{
		/* series of tests to determine which type of plot we will do */
		if( strcmp(plotCom.chPType[np]," MAP") == 0 )
		{
			/* thermal map */
			pltmap(np,chCall);
		}
		else if( strcmp(plotCom.chPType[np] ,"CONT") == 0 || 
			      strcmp(plotCom.chPType[np] ,"CRAW") == 0 || 
					strcmp(plotCom.chPType[np] ,"DIFF") == 0 || 
					strcmp(plotCom.chPType[np] ,"REFL") == 0 || 
					strcmp(plotCom.chPType[np] ,"EMIT") == 0 || 
					strcmp(plotCom.chPType[np] ,"CPHT") == 0 || 
					strcmp(plotCom.chPType[np] ,"OUTW") == 0 )
		{
			/* this is a contiuum plot of some kind */
			pltcon(np,chCall);
		}

		else if( 
			strcmp(plotCom.chPType[np] ,"OPAA") == 0 || 
			strcmp(plotCom.chPType[np] ,"OPAS") == 0 || 
			strcmp(plotCom.chPType[np] ,"OPAT") == 0 )
		{
			/*  absorption, scattering, or total opacity */
			pltopc(np,chCall);
		}
		else
		{
			fprintf( ioQQQ, " PLOT type=%4.4s not known.  STOP\n", 
			  plotCom.chPType[np] );
			cdEXIT(EXIT_FAILURE);
		}
	}

	return;
}

/*pltcon generate plot of continuum array */

STATIC void pltcon(
  long int np, 
  const char *chCall)
{
	char chSymPlt2[3], 
	  chXtitle[23];
	char chSym, 
	  chSymPlt1;
	long int i;
	double contin, 
	  ymin2;
	static double xmax, 
	  xmin, 
	  ymax, 
	  ymin;
	static realnum *y/*[rfield.nupper]*/, 
	  *y2/*[rfield.nupper]*/; 

	DEBUG_ENTRY( "pltcon()" );

	if( strcmp(chCall,"FIRST") == 0 )
	{ 
		return;
	}

	xmin = rfield.anulog[0];
	xmin = MAX2((double)plotCom.pltxmn[np],xmin);
	xmax = rfield.anulog[rfield.nflux-1];
	xmax = MIN2(xmax,(double)plotCom.pltxmx[np]);

	if( plotCom.lgPltTrace[np] )
	{
		fprintf( ioQQQ, " XMIN, XMAX=%12.4e%12.4e NFLUX=%4ld\n", 
		  xmin, xmax, rfield.nflux );
	}

	if( strcmp(plotCom.chPType[np],"REFL") == 0 && geometry.lgSphere )
	{
		fprintf( ioQQQ, " Reflected continuum not computed when SPHERE set.\n" );
		return;
	}

	y = (realnum*)MALLOC((size_t)rfield.nupper*sizeof(realnum) );
	y2 = (realnum*)MALLOC((size_t)rfield.nupper*sizeof(realnum) );

	/* these will be the default symbols for first and second plot */
	chSymPlt1 = '.';
	strcpy( chSymPlt2, "o " );
	ymin = FLT_MAX;
	ymin2 = FLT_MAX;
	ymax = -FLT_MAX;
	for( i=0; i < rfield.nflux; i++ )
	{
		if( (double)rfield.anulog[i] > xmin && (double)rfield.anulog[i] < xmax )
		{
			if( strcmp(plotCom.chPType[np],"CONT") == 0 )
			{
				y[i] = (realnum)log10(MAX2(rfield.flux_total_incident[0][i]/rfield.widflx[i]*
				  rfield.anu2[i],1e-37));
				/* >>chng 01 jul 13, add rfield.ConEmitReflec[0][i] */
				contin = rfield.flux[0][i] + rfield.ConEmitOut[0][i]*geometry.covgeo + rfield.ConEmitReflec[0][i];
				y2[i] = (realnum)MAX2((contin/rfield.widflx[i]+(rfield.outlin[0][i]+rfield.outlin_noplot[i])/
					rfield.anu[i]*geometry.covgeo)*
					rfield.anu2[i]*radius.r1r0sq,1e-37);
				y2[i] = (realnum)log10(y2[i]);
			}
			else if( strcmp(plotCom.chPType[np],"CPHT") == 0 )
			{
				/*    plot continuum as photons */
				y[i] = (realnum)log10(MAX2(rfield.flux_total_incident[0][i]/rfield.widflx[i],
				  1e-37));
				contin = rfield.flux[0][i] + rfield.ConEmitOut[0][i]*geometry.covgeo/rfield.widflx[i];
				y2[i] = (realnum)MAX2((contin+(rfield.outlin[0][i]+rfield.outlin[0][i])/rfield.anu[i]*
					geometry.covgeo)* radius.r1r0sq,1e-37);
				y2[i] = (realnum)log10(y2[i]);
			}
			else if( strcmp(plotCom.chPType[np],"REFL") == 0 )
			{
				/*    plot "reflected" continuum from last zone only */
				y[i] = (realnum)log10(MAX2((rfield.ConEmitReflec[0][i]/rfield.widflx[i]+
				  rfield.reflin[0][i])*rfield.anu2[i],1e-37));
				y2[i] = y[i];
			}
			else if( strcmp(plotCom.chPType[np],"EMIT") == 0 )
			{
				/*    plot "emitted" continuum from both sides of cloud */
				y[i] = (realnum)log10(MAX2(
					((rfield.ConEmitReflec[0][i]+rfield.ConEmitOut[0][i])/
					rfield.widflx[i]+
					(rfield.outlin[0][i]+rfield.outlin_noplot[i]+rfield.reflin[0][i])/rfield.anu[i] )*
					rfield.anu2[i],1e-37));
				y2[i] = y[i];
			}
			else if( strcmp(plotCom.chPType[np],"OUTW") == 0 )
			{
				/*    plot outward and attenuated incident continuum */
				chSymPlt1 = 'i';
				y[i] = (realnum)log10(MAX2(rfield.flux[0][i]*opac.opacity_abs[i],
				  1e-37));
				strcpy( chSymPlt2, "o " );
				y2[i] = (realnum)log10(MAX2((rfield.outlin[0][i]+rfield.outlin_noplot[i]+rfield.ConEmitOut[0][i])*
					opac.opacity_abs[i],1e-37));
			}
			else if( strcmp(plotCom.chPType[np],"DIFF") == 0 )
			{
				/*    plot "diffuse" continuum from last zone only */
				y[i] = (realnum)log10(MAX2(rfield.ConEmitLocal[nzone][i]*rfield.anu2[i]/
				  rfield.widflx[i],1e-37));
				y2[i] = y[i];
			}
			else if( strcmp(plotCom.chPType[np],"CRAW") == 0 )
			{
				y[i] = (realnum)log10(MAX2(rfield.flux_total_incident[0][i],1e-37));
				y2[i] = (realnum)MAX2((rfield.flux[0][i]+
				  rfield.otscon[i]+rfield.otslin[i]+rfield.outlin[0][i]+rfield.outlin_noplot[i]+
				  rfield.ConEmitOut[0][i])*radius.r1r0sq,1e-37);
				y2[i] = (realnum)log10(y2[i]);
			}

			if( y[i] > -36.9 )
			{
				ymin = MIN2(ymin,(double)y[i]);
			}

			if( y2[i] > -36.9 )
			{
				ymin2 = MIN2(ymin2,(double)y2[i]);
			}

			ymax = MAX2(ymax,(double)y[i]);
			ymax = MAX2(ymax,(double)y2[i]);
		}
	}

	if( trace.lgTrace )
	{
		fprintf( ioQQQ, " PLOT called for the first time, YMAX, MIN=%10.2e%10.2e\n", 
		  ymax, ymin );
	}

	/* lower min by at most 5 dex below peak */
	ymin2 = MAX3(ymax-5.,-35.,ymin2);

	/* make sure there is room at the bottom */
	ymin = MIN3(ymin2,ymin,ymax-1.);

	/* emitted continuum is thermal, so goes to zero */
	if( strcmp(plotCom.chPType[np],"EMIT") == 0 )
	{
		ymin = MAX2(ymin,ymax-4.);
	}

	if( plotCom.lgPltTrace[np] )
	{
		fprintf( ioQQQ, " YMAX, MIN=%14.4e%14.4e Npnts=%4ld\n", 
		  ymax
		  , ymin, rfield.nflux );
	}
	strcpy( chXtitle, "Log(nu fnu) vs LOG(nu)" );

	chSym = chSymPlt1;

	pltr(rfield.anulog,y,rfield.nflux,xmin,xmax,ymin,ymax,chSym,chXtitle
	  ,1,plotCom.lgPltTrace[np]);

	chSym = chSymPlt2[0];

	pltr(rfield.anulog,y2,rfield.nflux,xmin,xmax,ymin,ymax,chSym,chXtitle
	  ,3,plotCom.lgPltTrace[np]);

	free( y );
	free( y2 );
	return;
}

/*pltmap generate plot of heating and cooling map */

STATIC void pltmap(
  long int np, 
  const char *chCall)
{
	char chXtitle[23];
	static bool lgTlkSav;
	char chSym;

	long int i;

	static double xmax, 
	  xmin, 
	  ymax, 
	  ymin;

	DEBUG_ENTRY( "pltmap()" );

	if( strcmp(chCall,"FIRST") == 0 )
	{ 
		return;
	}

	lgTlkSav = called.lgTalk;
	called.lgTalk = false;
	hcmap.lgMapBeingDone = true;
	hcmap.RangeMap[0] = (realnum)pow((realnum)10.f,plotCom.pltxmn[np]);
	hcmap.RangeMap[1] = (realnum)pow((realnum)10.f,plotCom.pltxmx[np]);
	map_do(ioQQQ, " map");
	called.lgTalk = lgTlkSav;

	for( i=0; i < hcmap.nmap; i++ )
	{
		hcmap.temap[i] = (realnum)log10(hcmap.temap[i]);
	}

	xmin = MIN2(hcmap.temap[0],hcmap.temap[hcmap.nmap-1]);
	xmin = MAX2((double)plotCom.pltxmn[np],xmin);
	xmax = MAX2(hcmap.temap[0],hcmap.temap[hcmap.nmap-1]);
	xmax = MIN2(xmax,(double)plotCom.pltxmx[np]);

	if( plotCom.lgPltTrace[np] )
	{
		fprintf( ioQQQ, " xmin, xmax=%12.4e%12.4e nmap=%4ld\n", 
		  xmin, xmax, hcmap.nmap );
	}

	ymin = FLT_MAX;
	ymax = -FLT_MAX;

	for( i=0; i < hcmap.nmap; i++ )
	{
		if( (double)hcmap.temap[i] > xmin && (double)hcmap.temap[i] < xmax )
		{
			hcmap.hmap[i] = (realnum)log10(MAX2(hcmap.hmap[i],1e-35));
			hcmap.cmap[i] = (realnum)log10(MAX2(hcmap.cmap[i],1e-35));
			if( hcmap.cmap[i] > -34. )
			{
				ymin = MIN3(ymin,hcmap.hmap[i],hcmap.cmap[i]);
			}
			else
			{
				ymin = MIN2(ymin,(double)hcmap.hmap[i]);
			}
			ymax = MAX3(ymax,hcmap.hmap[i],hcmap.cmap[i]);
		}
	}

	if( trace.lgTrace )
	{
		fprintf( ioQQQ, " PLOT called for the first time, YMAX, MIN=%10.2e%10.2e\n", 
		  ymax, ymin );
	}

	if( plotCom.lgPltTrace[np] )
	{
		fprintf( ioQQQ, " YMAX, MIN=%14.4e%14.4e Npnts=%4ld\n", 
		  ymax, ymin, hcmap.nmap );
	}

	chSym = 'H';
	strcpy( chXtitle, "heating - cooling v te" );

	pltr(hcmap.temap,hcmap.hmap,hcmap.nmap,xmin,xmax,ymin,ymax,chSym,
	  chXtitle,1,plotCom.lgPltTrace[np]);

	chSym = 'C';

	pltr(hcmap.temap,hcmap.cmap,hcmap.nmap,xmin,xmax,ymin,ymax,chSym,
	  chXtitle,3,plotCom.lgPltTrace[np]);

	return;
}

/*pltopc generate plot of local gas opacity */
STATIC void pltopc(
  long int np, 
  const char *chCall)
{
	char chXtitle[23];
	char chSym;
	long int i;
	double arg1, 
	  arg2;
	static double xmax, 
	  xmin, 
	  ymax, 
	  ymin;
	static realnum *y/*[rfield.nupper]*/, 
	  *y2/*[rfield.nupper]*/;

	DEBUG_ENTRY( "pltopc()" );

	if( strcmp(chCall,"FIRST") == 0 )
	{ 
		return;
	}

	y = (realnum*)MALLOC((size_t)rfield.nupper*sizeof(realnum) );
	y2 = (realnum*)MALLOC((size_t)rfield.nupper*sizeof(realnum) );

	xmin = rfield.anulog[0];
	xmin = MAX2((double)plotCom.pltxmn[np],xmin);
	xmax = rfield.anulog[rfield.nflux-1];
	xmax = MIN2(xmax,(double)plotCom.pltxmx[np]);

	if( plotCom.lgPltTrace[np] )
	{
		fprintf( ioQQQ, " XMIN, XMAX=%12.4e%12.4e NFLUX=%4ld\n", 
		  xmin, xmax, rfield.nflux );
	}

	ymin = FLT_MAX;
	ymax = -FLT_MAX;

	for( i=0; i < rfield.nflux; i++ )
	{
		if( strcmp(plotCom.chPType[np],"OPAA") == 0 )
		{
			/*  absorption opacity */
			arg1 = opac.opacity_abs_savzon1[i];
			arg2 = opac.opacity_abs[i];
		}

		else if( strcmp(plotCom.chPType[np],"OPAS") == 0 )
		{
			/*  scattering opacity */
			arg1 = opac.opacity_sct_savzon1[i];
			arg2 = opac.opacity_sct[i];
		}

		else if( strcmp(plotCom.chPType[np],"OPAT") == 0 )
		{
			/*  total opacity */
			arg1 = opac.opacity_abs_savzon1[i] + opac.opacity_sct_savzon1[i];
			arg2 = opac.opacity_abs[i] + opac.opacity_sct[i];
		}

		else
		{
			/* this cannot happen since type was set to one of above */
			fprintf( ioQQQ, " pltopc type=%4.4s not known.  STOP\n", 
			  plotCom.chPType[np] );
			cdEXIT(EXIT_FAILURE);
		}

		y[i] = (realnum)log10(MAX2(arg1/dense.gas_phase[ipHYDROGEN],1e-35));
		y2[i] = (realnum)log10(MAX2(arg2/dense.gas_phase[ipHYDROGEN],1e-35));

		if( (double)rfield.anulog[i] > xmin && (double)rfield.anulog[i] < xmax )
		{
			ymin = MIN3(ymin,y[i],y2[i]);
			ymax = MAX3(ymax,y[i],y2[i]);
		}
	}

	if( trace.lgTrace )
	{
		fprintf( ioQQQ, " PLOT called for the first time, YMAX, MIN=%10.2e%10.2e\n", 
		  ymax, ymin );
	}

	/* lower min by factor of 10 to show absorption in next plot */
	ymin = MAX2(ymin-1.,-35.);
	ymax += 1.;
	if( plotCom.lgPltTrace[np] )
	{
		fprintf( ioQQQ, " YMAX, MIN=%14.4e%14.4e Npnts=%4ld\n", 
		  ymax, ymin, rfield.nflux );
	}

	strcpy( chXtitle, "Log(opacity) vs log(n)" );

	chSym = '.';
	pltr(rfield.anulog,y,rfield.nflux,xmin,xmax,ymin,ymax,chSym,chXtitle
	  ,1,plotCom.lgPltTrace[np]);

	chSym = 'o';
	pltr(rfield.anulog,y2,rfield.nflux,xmin,xmax,ymin,ymax,chSym,chXtitle
	  ,3,plotCom.lgPltTrace[np]);

	free(y);
	free(y2);
	return;
}

/*pltr core plotting routine for generating line printer plots */

STATIC void pltr(
	/* the x-axis */
	realnum x[], 
	/* the y-axi */
  realnum y[], 
  /* number of points */
  long int npnts, 
  /* mins and maxs, log of min and max of x-axis */
  double xmin, 
  double xmax, 
  double ymin, 
  double ymax, 
  /* plot symbol */
  char chSym, 
  char *chXtitle, 
  long int itim, 
  bool lgTrace)
{
	static char chPage[59][122];

	long int i, 
	  ix, 
	  iy, 
	  j, 
	  nc;

	/* the max number of decades we can plot */
#	define	NDECAD	18

	static long int jpnt[NDECAD], 
	  lowx, 
	  lx;

	static double xdec, 
	  xinc, 
	  ydown, 
	  yinc;

	/* this is log of smallestnumer in following set */
	const realnum xAxisMin = -8.f;

	static char chLab[NDECAD][5]={"1E-8","1E-7","1E-6","1E-5",
		"1E-4",".001","0.01"," 0.1","  1 ",
	  " 10 "," 100","1000","1E4 ","1E5 ","1E6 ","1E7 ","1E8 ","1E9 "};

	DEBUG_ENTRY( "pltr()" );

	/* ITIM=1, first call, =2 intermediate calls, =3 for last call*/
	if( itim == 1 )
	{
		/* first call, set left border of plot and clear out array */
		for( i=1; i < IHI; i++ )
		{
			chPage[i][0] = 'l';
			for( j=1; j < IWID; j++ )
			{
				chPage[i][j] = ' ';
			}
		}

		/* centered label for plot */
		strcpy( chPage[1], "                        " );
		strcat( chPage[1],  chXtitle );
		strcat( chPage[1], input.chTitle );

		/* one dex increments in x and y marked special */
		i = 1;
		ydown = 0.;
		yinc = (realnum)(IHI-2)/(ymax - ymin);
		nc = 0;

		while( i <= IHI && nc < 200 )
		{
			chPage[i-1][1] = '-';
			ydown += yinc;
			i = (long)(ydown + 1);
			nc += 1;
		}

		/* bottom increments of plot */
		for( i=0; i < IWID; i++ )
		{
			chPage[IHI-1][i] = '-';
		}

		if( xmin < xAxisMin )
		{
			fprintf(ioQQQ," plts: xmin is less than min value in array\n");
			cdEXIT(EXIT_FAILURE);
		}
		/* LX is pointer to label for number in x-axis in chLab */
		if( xmin < 0. )
		{
			lx = (long)(4.999-fabs(xmin));
			/* lx is the offset within the array of x-axis values */
			/* >>chng 99 jun 11 change to allow any min value of x-axis */
			lx = (long)(fabs(xAxisMin)-0.001-fabs(xmin));
			lx = MAX2(0,lx);
			/* this is lowest decade on plot */
			xdec = -floor(fabs(xmin)+1e-5);
		}
		else
		{
			double aa;
			lx = (long)MAX2(0.,4.+xmin);
			/* lx is the offset within the array of x-axis values */
			lx = (long)MAX2(0.,4.+xmin);
			/* >>chng 99 jun 11 change to allow any min value of x-axis */
			aa = fabs(xAxisMin);
			lx = (long)MAX2(0., aa-1. + xmin );
			xdec = floor(xmin+1e-5);
		}

		lowx = lx + 1;
		xinc = (realnum)(IWID-1)/(xmax - xmin);
		i = (long)MAX2(1.,(xdec-xmin)*xinc+1.);
		nc = 0;

		while( i < IWID && nc < 100 )
		{
			chPage[IHI-2][i - 1] = 'l';

			/*  fix position of x labels */
			lx = MIN2(lx+1,NDECAD);

			/*  slight offset to center label */
			jpnt[lx-1] = MAX2(0,i-3);
			jpnt[lx-1] = MIN2((long)IWID-4,jpnt[lx-1]);
			xdec += 1.;
			i = (long)MAX2(1.,(xdec-xmin)*xinc+1.);
			nc += 1;
		}
	}

	/* everything falls down through here */
	/* now fill in data, symbol is chSym */
	for( i=0; i < npnts; i++ )
	{
		if( (double)x[i] > xmin && (double)x[i] < xmax )
		{
			iy = (long)(IHI - MAX2(y[i]-ymin,0.)*yinc);
			iy = MAX2(1,iy);
			ix = (long)((x[i] - xmin)*xinc + 1);

			if( lgTrace )
			{
				fprintf( ioQQQ, " x, y, ix, iy=%7.3f%7.3f%4ld%4ld\n", 
				  x[i], y[i], ix, iy );
			}
			chPage[iy-1][ix - 1] = chSym;
		}
	}

	if( itim == 3 )
	{
		/* make the plot */
		fprintf( ioQQQ, "1\n" );
		for( i=1; i < IHI; i++ )
		{
			fprintf( ioQQQ, "     %121.121s\n", chPage[i] );
		}

		/* now put on label for X-axis */
		for( i=0; i < IWID; i++ )
		{
			chPage[0][i] = ' ';
		}

		for( i=lowx-1; i < lx; i++ )
		{
			/* copy the four char of the numeric string */
			strncpy(chPage[0]+jpnt[i] , chLab[i+1] , 4);
		}
		fprintf( ioQQQ, "     %121.121s\n", chPage[0] );
	}
	return;
}