Difference between revisions of "Plot skymap (xfig example)"
Line 1: | Line 1: | ||
[[File:Skymapcirclegr.png]] | [[File:Skymapcirclegr.png]] | ||
− | + | [code] | |
require("isisscripts"); | require("isisscripts"); | ||
Line 253: | Line 253: | ||
% Render | % Render | ||
p.render("source_skyplot.pdf"); | p.render("source_skyplot.pdf"); | ||
+ | [/code] |
Revision as of 15:59, 5 September 2018
[code] 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"); [/code]