/*plot master routine to generate some sort of plot */
#include <float.h>
#include <string.h>
#include "cddefines.h"
#include "itercnt.h"
#include "called.h"
#define	IHI	59
#define	IWID	121
#include "input.h"
#include "rfield.h"
#include "flxsav.h"
#include "trace.h"
#include "radius.h"
#include "sphere.h"
#include "sumoutinc.h"
#include "opac.h"
#include "mappar.h"
#include "punt.h"
#include "cdminmax.h"/* variable args */
#include "phycon.h"
#include "plot.h"
/*pltcon generate plot of continuum array */
void pltcon(long int np, 
  char *chCall);

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

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

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


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

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

	/* return if this is not the last iteration, or a plot not required,
	 * or we are not speaking */
	if( !plotCom.lgPlotON || !called.lgTalk )
	{ 
#		ifdef DEBUG_FUN
		fputs( " <->plot()\n", debug_fp );
#		endif
		return;
	}

	if( !IterCnt.lgLastIt && (strcmp(chCall,"FIRST") != 0) )
	{ 
#		ifdef DEBUG_FUN
		fputs( " <->plot()\n", debug_fp );
#		endif
		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] );
			puts( "[Stop in plot]" );
			exit(1);
		}
	}


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

/*pltcon generate plot of continuum array */

void pltcon(long int np, 
  char *chCall)
{
	char chSymPlt2[3], 
	  chXtitle[23];
	char chSym, 
	  chSymPlt1;
	long int i;
	double contin, 
	  ymin2;
	static double xmax, 
	  xmin, 
	  ymax, 
	  ymin;
	static float y[NCELL], 
	  y2[NCELL]; 

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


	if( strcmp(chCall,"FIRST") == 0 )
	{ 
#		ifdef DEBUG_FUN
		fputs( " <->pltcon()\n", debug_fp );
#		endif
		return;
	}

	xmin = rfield.anulog[0];
	xmin = MAX2(plotCom.pltxmn[np],xmin);
	xmax = rfield.anulog[rfield.nflux-1];
	xmax = MIN2(xmax,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 && sphere.lgSphere )
	{
		fprintf( ioQQQ, " Reflected continuum not computed when SPHERE set.\n" );
		
#		ifdef DEBUG_FUN
		fputs( " <->pltcon()\n", debug_fp );
#		endif
		return;
	}

	/* 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( rfield.anulog[i] > xmin && rfield.anulog[i] < xmax )
		{
			if( strcmp(plotCom.chPType[np],"CONT") == 0 )
			{
				y[i] = (float)log10(MAX2(flxsav.FluxSave[i]/rfield.widflx[i]*
				  rfield.anu2[i],1e-37));
				contin = rfield.flux[i] + (rfield.ConOutNoInter[i] + 
				  rfield.outcon[i])*sphere.covgeo;
				y2[i] = (float)MAX2((contin/rfield.widflx[i]+rfield.outlin[i]/
					rfield.anu[i]*sphere.covgeo)*
					rfield.anu2[i]*radius.r1r0sq,1e-37);
				y2[i] = (float)log10(y2[i]);
			}
			else if( strcmp(plotCom.chPType[np],"CPHT") == 0 )
			{
				/*    plot continuum as photons */
				y[i] = (float)log10(MAX2(flxsav.FluxSave[i]/rfield.widflx[i],
				  1e-37));
				contin = rfield.flux[i] + (rfield.ConOutNoInter[i] + 
				  rfield.outcon[i])*sphere.covgeo/rfield.widflx[i];
				y2[i] = (float)MAX2((contin+rfield.outlin[i]/rfield.anu[i]*
					sphere.covgeo)* radius.r1r0sq,1e-37);
				y2[i] = (float)log10(y2[i]);
			}
			else if( strcmp(plotCom.chPType[np],"REFL") == 0 )
			{
				/*    plot "reflected" continuum from last zone only */
				y[i] = (float)log10(MAX2((rfield.refcon[i]/rfield.widflx[i]+
				  rfield.reflin[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] = (float)log10(MAX2(
					((rfield.refcon[i]+rfield.ConOutNoInter[i]+rfield.outcon[i])/
					rfield.widflx[i]+
					(rfield.outlin[i]+rfield.reflin[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] = (float)log10(MAX2(rfield.flux[i]*opac.opac[i],
				  1e-37));
				strcpy( chSymPlt2, "o " );
				y2[i] = (float)log10(MAX2((rfield.outlin[i]+rfield.outcon[i]-
				  SumOutInc.SavOutCon[i])*opac.opac[i],1e-37));
			}
			else if( strcmp(plotCom.chPType[np],"DIFF") == 0 )
			{
				/*    plot "diffuse" continuum from last zone only */
				y[i] = (float)log10(MAX2(rfield.DiffLocal[i]*rfield.anu2[i]/
				  rfield.widflx[i],1e-37));
				y2[i] = y[i];
			}
			else if( strcmp(plotCom.chPType[np],"CRAW") == 0 )
			{
				y[i] = (float)log10(MAX2(flxsav.FluxSave[i],1e-37));
				y2[i] = (float)MAX2((rfield.flux[i]+rfield.ConOutNoInter[i]+
				  rfield.otscon[i]+rfield.otslin[i]+rfield.outlin[i]+
				  rfield.outcon[i])*radius.r1r0sq,1e-37);
				y2[i] = (float)log10(y2[i]);
			}

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

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

			ymax = MAX2(ymax,y[i]);
			ymax = MAX2(ymax,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 = vfmax(ymax-5.,-35.,ymin2,FEND);

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

	/* 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]);

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

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

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

	long int i;

	static double xmax, 
	  xmin, 
	  ymax, 
	  ymin;

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

	if( strcmp(chCall,"FIRST") == 0 )
	{ 
#		ifdef DEBUG_FUN
		fputs( " <->pltmap()\n", debug_fp );
#		endif
		return;
	}

	lgTlkSav = called.lgTalk;
	called.lgTalk = FALSE;
	MapPar.RangeMap[0] = (float)pow(10.,plotCom.pltxmn[np]);
	MapPar.RangeMap[1] = (float)pow(10.,plotCom.pltxmx[np]);
	punt(" map");
	called.lgTalk = lgTlkSav;

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

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

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

	ymin = FLT_MAX;
	ymax = -FLT_MAX;

	for( i=0; i < MapPar.nmap; i++ )
	{
		if( MapPar.temap[i] > xmin && MapPar.temap[i] < xmax )
		{
			MapPar.hmap[i] = (float)log10(MAX2(MapPar.hmap[i],1e-35));
			MapPar.cmap[i] = (float)log10(MAX2(MapPar.cmap[i],1e-35));
			if( MapPar.cmap[i] > -34. )
			{
				ymin = vfmin(ymin,MapPar.hmap[i],MapPar.cmap[i],FEND);
			}
			else
			{
				ymin = MIN2(ymin,MapPar.hmap[i]);
			}
			ymax = vfmax(ymax,MapPar.hmap[i],MapPar.cmap[i],FEND);
		}
	}

	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, MapPar.nmap );
	}

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

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

	chSym = 'C';

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


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

/*pltopc generate plot of local gas opacity */

void pltopc(long int np, 
  char *chCall)
{
	char chXtitle[23];
	char chSym;
	long int i;
	double arg1, 
	  arg2;
	static double xmax, 
	  xmin, 
	  ymax, 
	  ymin;
	static float y[NCELL], 
	  y2[NCELL];

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

	if( strcmp(chCall,"FIRST") == 0 )
		{ 
#	ifdef DEBUG_FUN
		fputs( " <->pltopc()\n", debug_fp );
#	endif
		return;}
	xmin = rfield.anulog[0];
	xmin = MAX2(plotCom.pltxmn[np],xmin);
	xmax = rfield.anulog[rfield.nflux-1];
	xmax = MIN2(xmax,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.opacsv[i];
			arg2 = opac.opac[i];
		}

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

		else if( strcmp(plotCom.chPType[np],"OPAT") == 0 )
		{
			/*  total opacity */
			arg1 = opac.opacsv[i] + opac.sctsav[i];
			arg2 = opac.opac[i] + opac.scatop[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] );
			puts( "[Stop in pltopc]" );
			exit(1);
		}

		y[i] = (float)log10(MAX2(arg1/phycon.hden,1e-35));
		y2[i] = (float)log10(MAX2(arg2/phycon.hden,1e-35));

		if( rfield.anulog[i] > xmin && rfield.anulog[i] < xmax )
		{
			ymin = vfmin(ymin,y[i],y2[i],FEND);
			ymax = vfmax(ymax,y[i],y2[i],FEND);
		}
	}

	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]);


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

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

void pltr(
	/* the x-axis */
	float x[], 
	/* the y-axi */
  float 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 */
  long int chSym, 
  char *chXtitle, 
  long int itim, 
  int 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 float 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 "};

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

	/* 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 = (float)(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");
			puts( "[Stop in ParseTable]" );
			exit(1);
		}
		/* 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 = (long)MAX2(0.,lx);
			/* this is lowest decade on plot */
			xdec = -(long)(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 = (long)(xmin+1e-5);
		}

		lowx = lx + 1;
		xinc = (float)(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(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( x[i] > xmin && 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] = (char)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] );
	}

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