Continuum functions (tikz plot)

S-Lang code:

#!/usr/bin/env isis

require("isisscripts");
require("tikz");
use_localmodel("pulsar");
tikz_add_default_latex_package("txfonts");


% functions to check: cutoffpl, highecut, fdcut

variable funcs=["cutoffpl(1)","powerlaw(1)*highecut(1)","powerlaw(1)*fdcut(1)"];
% PHO FOLD (and only if) CUT
variable varPars={[2,3],[2,4,3],[2,4,3]};;
variable parNames=["\\Gamma","E_\\mathrm{fold}","E_\\mathrm{cutoff}"];
variable unit=["","keV","keV"];
variable lower=[-1,1,10];
variable upper=[1,6,20];

variable whereFOLD=[3,4,4];
variable whereCUT=[-1,3,3];


variable funcnames=["\\texttt{cutoffpl}","\\texttt{highecut}","\\texttt{FDcut}"];
variable formulas=["","",""];
formulas[0]="$K \\cdot E^{-\\Gamma} \\exp\\left(-\\frac{E}{E_\\mathrm{fold}}\\right)$";
formulas[1]="$K \\cdot E^{-\\Gamma} \\exp\\left(-\\frac{E-E_\\mathrm{cutoff}}{E_\\mathrm{fold}}\\right)$";
formulas[2]="$K \\cdot E^{-\\Gamma}\\left[ 1 + \\exp \\left( \\frac{E-E_\\mathrm{cutoff}}{E_\\mathrm{fold}}\\right) \\right]^{-1}$";

variable numSteps=8;
variable func,varPar,ii,jj,kk,lo,hi;
variable Emin=1;
variable Emax=80;

%energy grid
(lo,hi)=log_grid(Emin,Emax,1000);
variable middle=0.5*(lo+hi);
variable width=hi-lo;


variable pho=-0.5;
variable ecut=15;
variable efold=5;

variable compPlot=Struct_Type[length(funcs)];

for(ii=0;ii<length(funcs);ii++){
	func=funcs[ii];
	fit_fun(func);


	% create structure that holds this column (length is always free since
	% cutoffpl gets the info tab
	variable column=Struct_Type[3];
	%loop over variable params
	for(jj=0;jj<length(varPars[ii]);jj++){
		
		% establish basis params
		set_par(1,1,1); %norm set to 1
		if(ii!=0){set_par(whereCUT[ii],ecut);};  
		set_par(whereFOLD[ii],efold);
		set_par(2,pho);
		% basic setup complete! now vary the rest
	
		% add a subplot
		column[jj]=tikz_plot_new(5,5);
		column[jj].world(Emin,Emax,5e-2,2e1;xlog,ylog);

		% decativate x1axis, activate x2axis
		column[jj].x1axis(;format="");
		column[jj].x2axis(;on);
		
		% decativate y1axis, activate y2axis
		column[jj].y1axis(;format="");
		%activate only for last column
		if(ii==2){
			column[jj].y2axis(;on,format="$10^{%.0f}$");
			column[jj].y2label("Flux density [arb. units]");
		};

		% slope of parameter change	
		variable slope=1.0*(upper[jj]-lower[jj])/(numSteps-1);
	
		% loop over all the different values to try!
		for(kk=0;kk<numSteps;kk++){

			% determine the parameter
			set_par(varPars[ii][jj],slope*kk+lower[jj]);
			variable fluxDen=eval_fun_keV(lo,hi)/width;

			% determine the plot color
			tikz_new_color("plotcol",sprintf("blue!%d!red",100*kk/numSteps));
		
			% plot the obtained curve
			column[jj].plot(middle,fluxDen;color="plotcol");
		};


	if((ii==0 && jj<2) || ii!=0){
		column[jj].xylabel(0.04,0.9,sprintf("$%s={\\color{red}%d}$--${\\color{blue}%d}$\\,%s",parNames[jj],lower[jj],upper[jj],unit[jj]);world0,anchor="west");
	};


	% add dashed line for the cutoff energy
	if(ii>0 && jj<2){
		column[jj].plot([ecut,ecut],[0,1];world10,line=1);
	}



	};

	% add the info block and parameter border info	
	if(ii==0){
		column[2]=tikz_plot_new(5,5);
		column[2].xaxis(;off);
		column[1].xaxis(;on);
		column[2].yaxis(;off);
		column[2].xylabel(0,0.8,"\\textbf{Constant parameters}";anchor="west",world0);
		variable tablehead="\\begin{tabular}{lr} \\hline Parameter & Value \\\\ \\hline";
		variable table=sprintf("%s Normalization & 1 \\\\ Photon Index $\\Gamma$& $%.1f$ \\\\ $E_\\mathrm{cutoff}$ & $%.0f$\\,keV \\\\ $E_\\mathrm{fold}$ & $%.0f$\\,keV \\\\ \\hline  \\end{tabular}",tablehead,pho,ecut,efold);
		column[2].xylabel(0,0.5,table;anchor="west",world0);
	};
	
	compPlot[ii]=tikz_multiplot(__push_array(column);x2label="Energy [keV]",title=sprintf("\\begin{tabular}{c} %s \\\\ %s \\end{tabular}",funcnames[ii],formulas[ii]));	
};

% put it all together and render
variable prep=tikz_new_hbox_compound(__push_array(compPlot),-0.1;just=1);



prep.render("continuum_functions.pdf");

File:Continuum functions.pdf