Plot skymap (xfig example)

require("isisscripts");

%runs in /home/krauss/scripts/neutrinos100tev
% see there for the files used below

%%%%%%%%%%%%%%%%%%%%%%%%%
% Read coordiates
variable abc = ascii_read_table("nuev.txt",[{"%s","id"},{"%F","e"},{"%F","ra"},{"%F","dec"},{"%F","angr"}]);


%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot setup

xfig_new_color ("galactic_col", 0x0000FF);
xfig_new_color ("gridcol",      0x777777);
variable N = 300; % #steps for axes
variable plot_width  = 16;
variable plot_height = 8;
variable p = xfig_plot_new(plot_width,plot_height);
p.world (-2,2,-1,1);
p.axis(;off);



%%%%%%%%%%%%%%%%%%%%%%%%%
% GRID

variable RA, DEC;
variable ra_adjust  = 0.6;
variable dec_adjust = 0.6;

foreach RA ([-180:180:60])
{
  p.plot(Hammer_projection ( Double_Type[N]+RA , [-90:90:#N]; deg, normalized) ; color="gridcol", width=2);
  p.xylabel(Hammer_projection ( RA, 0 ; deg, normalized),sprintf("%dh",((360-RA)/15) mod 24), ra_adjust, dec_adjust);
}
ra_adjust  = -0.6;
dec_adjust = 0.6;
foreach DEC ([-60:60:30])
{
  p.plot(Hammer_projection ( [-180:180:#N], Double_Type[N]+DEC ; deg, normalized) ; color="gridcol", width=2);
  p.xylabel(Hammer_projection ( 0,  DEC ; deg, normalized),sprintf("%d$^\\circ$",DEC), ra_adjust, dec_adjust);
}

%%%%%%%%%%%%%%%%%%%%%%%%%
% EVENTS crosses

p.plot(Hammer_projection (  -abc.ra,  abc.dec ; deg, normalized); sym = "+", fill = 20, size = 0.7, width=3);

variable i;
_for i(0, length(abc.dec)-1,1)
{
    p.xylabel(Hammer_projection (  -abc.ra[i],  abc.dec[i] ; deg, normalized), "\scriptsize "R+abc.id[i],-0.2,-0.8;depth=1);
}


p.plot(0.01,0.94; world0, sym = "+", fill = 20, size = 0.7);
p.xylabel(0.03,0.94,"IC events",-0.5,0; world0);

variable R=0.01;
variable tmp = [0:2*PI:#100];

p.plot((R*cos(tmp))+0.01, R*sin(tmp)+0.86;world00,line=0, color="black" );
p.xylabel(0.03,0.86,"R$_{50}$"R,-0.5,0; world0);


%%%%%%%%%%%%%%%%%%%%%%%%%
% Galactic plane

variable ra,dec,i;
variable galactic_col = "gal_plane";
(ra,dec) = RAdec_from_galLB ([0:360:#300],Double_Type[300]);
i = array_sort(pos_modulo(ra-180,360));
p.plot(Hammer_projection (-ra[0],dec[0]; deg, normalized); color="galactic_col", depth=100, width=2, sym = "x");
ra = ra[i]; dec = dec[i];
p.plot(Hammer_projection (-ra,dec; deg, normalized); color="galactic_col", depth=10, width=3);
p.xylabel(0.99,0.94,"Galactic Plane",0.5,0; world0,color =  "galactic_col");


p.plot(0.85,0.06; world0, sym = "x", color="gray", size = 0.7);
p.xylabel(0.87,0.06,"3FGL sources",-0.5,0; world0, color="gray");



%%%%%%%%%%%%%%%%%%%%%%%%%
% Circles

variable R; %angular resolution
define get_circle (R,ra,dec)
{
    ra = ra*PI/180;
    dec = dec*PI/180+PI/2;
    R = R* PI/180;

    %print(ra);
    %print(dec);
    variable tmp = [0:2*PI:#100];
    
    variable x = (sin(R)*cos(ra)*cos(dec)*cos(tmp)) -
    (sin(R)*sin(ra)*sin(tmp)) + (cos(R)*sin(dec)*cos(ra));
    
    variable y = (sin(R)*sin(ra)*cos(dec)*cos(tmp)) +
    (sin(R)*cos(ra)*sin(tmp)) + (cos(R)*sin(dec)*sin(ra));

    variable z = -sin(R)*sin(dec)*cos(tmp)+cos(R)*cos(dec);

    variable decnew = acos(z);
    variable ranew = atan2(y,x);
    %print(max(ranew));
    %print(max(decnew));
    return ranew*180/PI, (decnew-PI/2)*180/PI;
}


variable j;

_for i (0, length(abc.angr)-1,1)
{
    
    R=abc.angr[i];
    variable ranew,decnew;
    (ranew,decnew) = get_circle (R,abc.ra[i],abc.dec[i]);
    
    %print(R);
    
    variable indexbreak=0;
    variable breakflag=0;
    variable prevval = Double_Type[0];
    variable breakval = Double_Type[0];
    

    if (i==1)
    {
	%print(ranew);
    }
    _for j (0, length(ranew)-2,1)
    {


	
	if (abs(ranew[j]-ranew[j+1]) > 10)
	{

	    prevval = [prevval,ranew[j]];
	    breakval = [breakval,ranew[j+1]];
	    indexbreak=[indexbreak,j+1];
	    breakflag=1;
	}


    }


    if (breakflag == 1)
    {


	variable xnew = struct {ra = ranew, dec=decnew };
	
	variable xnew1 = struct_filter(xnew,where(xnew.ra <= 0);copy);
	variable xnew2 = struct_filter(xnew,where(xnew.ra >= 0);copy);


	
	if (abc.e[i]>1000)
	{
	    p.plot(Hammer_projection(-(xnew1.ra),xnew1.dec;deg, normalized) ;
	    line=0, color="gray" );
	    p.plot(Hammer_projection(-(xnew2.ra),xnew2.dec;deg, normalized) ;
	    line=0, color="gray" );
	}
	else
	{
	    p.plot(Hammer_projection(-xnew1.ra,xnew1.dec;deg, normalized) ;
	    line=0, color="red" );
	    
	    p.plot(Hammer_projection(-xnew2.ra,xnew2.dec;deg, normalized) ;
	    line=0, color="red" );
	}
    }
    else
    {
	if (abc.e[i]>1000)
	{
	    p.plot(Hammer_projection(-(ranew),decnew;deg, normalized) ;
	    line=0, color="gray" );
	}
	else
	{
	    p.plot(Hammer_projection(-ranew,decnew;deg, normalized) ;
	    line=0, color="red" );
	}
}

} 





%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 3FGL SRC
require("/home/krauss/scripts/neutrinos100tev/read3fgl.sl");
variable effpath = "/home/krauss/important_data/neutrinos_2015_12/";

variable event = [2,4,12,13,17,22,26,30,33,38,39,40,45,46,48,52];
variable j;
_for i (0, length(event)-1,1)
{

    variable b = ascii_read_table (effpath+"eventsummary2.txt",[{"%u","id"},{"%f","energy"},{"%s","b1"},{"%f","energyerr1"},{"%f","energyerr2"},{"%s","b2"},{"%s","time"},{"%f","dec"},{"%f","ra"},{"%s","ang"},{"%s","top"}]);
    
    struct_filter(b, where(b.id == event[i]));
    
    if (length(b.ra) > 1 || length(b.ra) == 0)
    {
	vmessage("Something went wrong, no or multiple events found!");
    }

    if (b.ang[0] == "<1.4")
    {
	variable ang = 1.4;
    }
    else
    {
	ang = atof(b.ang[0]);
    }
    
    variable all_src, evt;
    
    (evt, all_src) = get_fgl (b.id[0], b.ra[0], b.dec[0], ang);

    _for j(0, length(all_src.src)-1,1)
    {
	
	variable ratem = strchop(all_src.ra[j],' ',0);
	variable dectem = strchop(all_src.dec[j],' ',0);

	variable ranew = hr2ra(atof(ratem[1]),atof(ratem[2]),atof(ratem[3]));
	variable decnew = hr2dec(atof(dectem[1]),atof(dectem[2]),atof(dectem[3]));


	p.plot(Hammer_projection(-ranew,decnew;deg,normalized);color="gray",sym="x",size=0.35);

    }

}

%%%%%%%%%%%%%
% Render
p.render("source_skyplot.pdf");