Commit cbb77131 authored by Michal Bursa's avatar Michal Bursa
Browse files

added plots

parent 0ad3dc35
/data
/tlusty/*
!/tlusty/t205*.f
!/tlusty/data
!/tlusty/*.FOR
\ No newline at end of file
!/tlusty/*.FOR
/plots/*
!/plots/*.asy
!/plots/*.py
import graph;
import palette;
bounds shade_image(picture pic=currentpicture, pair[] z, real[] f,
range range=Full, pen[] palette)
{
if(z.length != f.length)
abort("z and f arrays have different lengths");
real m=min(f);
real M=max(f);
bounds bounds=range(pic,m,M);
real rmin=pic.scale.z.T(bounds.min);
real rmax=pic.scale.z.T(bounds.max);
palette=adjust(pic,m,M,rmin,rmax,palette);
// Crop data to allowed range and scale
if(range != Full || pic.scale.z.scale.T != identity ||
pic.scale.z.postscale.T != identity) {
scalefcn T=pic.scale.z.T;
real m=bounds.min;
real M=bounds.max;
f=map(new real(real x) {return T(min(max(x,m),M));},f);
}
int[] edges={0,0,1};
int N=palette.length-1;
int[][] trn=triangulate(z);
real step=rmax == rmin ? 0.0 : N/(rmax-rmin);
for(int i=0; i < trn.length; ++i) {
int[] trni=trn[i];
int i0=trni[0], i1=trni[1], i2=trni[2];
pen color(int i) {return palette[round((f[i]-rmin)*step)];}
gouraudshade(pic,z[i0]--z[i1]--z[i2]--cycle,
new pen[] {color(i0),color(i1),color(i2)},edges);
}
return bounds; // Return bounds used for color space
}
string filename = "../a000/l010/0050/stage1.15i_ang";
file f = input(filename).line();
real[][] data = f.dimension(0,0);
close(f);
pair[] xy;
real[] f;
for (int i=0; i<data.length; ++i) {
real theta = data[i][0]/180*pi; // read first value (value of theta angle)
data[i].delete(0); // remove first value from array
pair[] theta_xy = sequence(
new pair(int j) {
real phi = (j+0.5) * 2*pi/data[i].length;
real r = theta/(pi/2);
return (500*r*cos(phi), 500*r*sin(phi));
},
data[i].length
);
real[] theta_f = sequence(
new real(int j) {
return (data[i][j]<=0.0) ? data[i][j] : log10(data[i][j]);
},
data[i].length
);
xy.append(theta_xy);
f.append(theta_f);
}
//write(data.length);
//write(xy.length);
write(xy.length);
while (xy.length > 1000) {
int[][] trn=triangulate(xy);
// step 1 - find triangle with least difference
int[] toremove;
for (int i=0; i<trn.length; ++i) {
int i0=trn[i][0], i1=trn[i][1], i2=trn[i][2];
real diff = max(f[i0], f[i1], f[i2]);
if (diff == 0.0) {
pair p12 = 0.5(xy[i1]+xy[i2]);
pair t = 1/3*xy[i0] + 2/3*p12;
real val = f[i0];
//xy.push(t); f.push(val);
toremove.push(i0);
toremove.push(i1);
toremove.push(i2);
}
}
write(xy.length);
toremove = sort(toremove);
int lastremoved = -1;
for (int i=toremove.length-1; i>=0; --i) {
if (lastremoved == toremove[i]) continue;
xy.delete(toremove[i]);
f.delete(toremove[i]);
lastremoved = toremove[i];
}
write(xy.length);
//abort("iter 0");
//continue;
break;
// step 1 - find triangle with least difference
real mindiff = 1e200;
int imindiff = -1;
int[][] trn=triangulate(xy);
for(int i=0; i<trn.length; ++i) {
int i0=trn[i][0], i1=trn[i][1], i2=trn[i][2];
real diff = max(f[i0], f[i1], f[i2]);
if (diff < mindiff) {
mindiff = diff;
imindiff = i;
}
}
// step 2 - remove one point from the array
/*
int[] trni=trn[imindiff];
int i0=trni[0], i1=trni[1], i2=trni[2];
int iremove;
if ((length(xy[i0])>=length(xy[i1])) && (length(xy[i0])>=length(xy[i2]))) {
iremove = (abs(f[i0]-f[i1]) < abs(f[i0]-f[i2])) ? 1 : 2;
} else
if ((length(xy[i1])>=length(xy[i0])) && (length(xy[i1])>=length(xy[i2]))) {
iremove = (abs(f[i1]-f[i0]) < abs(f[i1]-f[i2])) ? 0 : 2;
} else
if ((length(xy[i2])>=length(xy[i0])) && (length(xy[i2])>=length(xy[i1]))) {
iremove = (abs(f[i2]-f[i0]) < abs(f[i2]-f[i1])) ? 0 : 1;
}
*/
int[] trni=trn[imindiff];
int i0=trni[0], i1=trni[1], i2=trni[2];
write(imindiff, mindiff);
int iremove;
real d01 = length(xy[i0]-xy[i1]);
real d02 = length(xy[i0]-xy[i2]);
real d12 = length(xy[i1]-xy[i2]);
if ((d01<=d02) && (d01<=d12)) {
iremove = (d02<d12) ? 0 : 1;
} else
if ((d02<=d12) && (d02<=d01)) {
iremove = (d01<d12) ? 0 : 2;
} else
if ((d12<=d01) && (d12<=d02)) {
iremove = (d01<d02) ? 1 : 2;
}
xy.delete(iremove);
f.delete(iremove);
write(xy.length);
}
//write(xy.length);
bounds range = shade_image(xy, f, Full, Rainbow());
//palette(bar,"$A$",range,(0,0),(0.5cm,8cm),Right,Palette,
// PaletteTicks("$%+#.1f$"));
import graph;
real frequency_integral(real[] E, real[] I) {
int N = E.length;
real flux = 0.0;
flux += I[0] * 0.5*abs(E[1]-E[0]);
for (int i=1; i<N-1; ++i) flux += I[i] * 0.5*abs(E[i+1]-E[i-1]);
flux += I[N-1] * 0.5*abs(E[N-2]-E[N-1]);
return flux;
}
void read_spectrum(string filename, real[] A, real[] I)
{
A.delete();
I.delete();
file f = input(filename).line();
real[][] raw = f.dimension(0,0);
close(f);
// if row[0] contains only one number, we have module 15 and we skip first 3 lines
// if row[0] contains two numbers, we have module 15i and we skip first 2 lines
int row_offset = (raw[0].length == 1) ? 3 : 2;
real[] cosines = raw[1];
real[][] tmp = transpose(raw[row_offset:]); // input radiances in [erg/cm2/s/Hz/srad]
for (int i=0; i<cosines.length; ++i) {
A[i] = acos(cosines[i]);
I[i] = log10(frequency_integral(tmp[0], tmp[2+i]));
if (I[i]<0.0) I[i]=0.0;
}
write(A*180/pi);
if (A[0] < A[1]) {
int e = A.length-1;
real t;
for(int i=0; i<A.length/2; ++i) {
t = A[i]; A[i] = A[e]; A[e] = t;
t = I[i]; I[i] = I[e]; I[e] = t;
e=e-1;
}
}
write(A*180/pi);
}
// --- settings ----------------
//real mass = 10;
//real spin = 0.65;
//real alpha = 0.100;
//real mdot[] = {0.1, 0.3, 0.6};
//pen colors[] = {black,deepblue,heavyred}; if (settings.gray) colors = new pen[] {black,black,black};
//real inc = 66;
real xmin = 0.01;
real xmax = 30.0;
real ymin = -infinity;
real ymax = infinity;
real size_factor = 1.0;
real size_aspect = 3/2;
defaultpen(currentpen+linewidth(1pt));
usersetting();
pen pen_original = black+linetype("3 3");
pen pen_retrad = black+linetype("0 3");
pen pen_final = black+solid;
//------ plot graph --------------
picture plot_spectrum() {
picture pic;
size(pic, 80mm*size_factor*size_aspect, 80mm*size_factor, Aspect);
scale(pic, Linear, Linear);
real[] A;
real[] I;
read_spectrum("../a000/l010/0130/stage1.15", A, I);
(A = pi/2-A).append(reverse(pi-A));
(I = I).append(reverse(I));
draw(pic, polargraph(pic, I, A, operator ..), pen_original, "original");
read_spectrum("../a000/l010/0130/stage1.15i", A, I);
(A = pi/2-A).append(reverse(pi-A));
(I = I).append(reverse(I));
draw(pic, polargraph(pic, I, A, operator --), pen_retrad, "original");
read_spectrum("../a000/l010/0130/stage2.15", A, I);
(A = pi/2-A).append(reverse(pi-A));
(I = I).append(reverse(I));
draw(pic, polargraph(pic, I, A, operator ..), pen_final, "original");
// read_spectrum("../a000/l010/0013/stage1.15i", E, I);
// draw(pic, graph(pic, E, I), pen_retrad, "retrad");
// read_spectrum("../a000/l010/0013/stage2.15", E, I);
// draw(pic, graph(pic, E, I), pen_final, "final");
// real ymax = 10^(round(log10(max(I)))+0.5);
// real ymin = ymax*1e-5;
// label flux curves
// label(pic, rotate(-62)*"$L=0.1$", (log10(9.5),log10(5*ymin)), W, Helvetica()+fontsize(9pt), NoFill);
// label(pic, rotate(-62)*"$L=0.3$", (log10(13.5),log10(5*ymin)), W, Helvetica()+fontsize(9pt), NoFill);
// label(pic, rotate(-62)*"$L=0.6$", (log10(21),log10(5*ymin)), W, Helvetica()+fontsize(9pt), NoFill);
//xlimits(pic, xmin, xmax);
//ylimits(pic, 0, infinity);
//crop(pic);
xaxis(pic, "energy [keV]", BottomTop, LeftTicks(DefaultFormat), above=true);
yaxis(pic, "flux [erg/cm$^2$/keV/s]", LeftRight, RightTicks, above=true);
/*
// attach legend
legendlinelength = 40;
legendmargin = 6;
Legend l1 = Legend(" NT", black+linewidth(1pt)+dashed+Helvetica()+fontsize(9pt));
Legend l2 = Legend(" SD", black+linewidth(1pt)+solid+Helvetica()+fontsize(9pt));
pic.legend.push(l1);
pic.legend.push(l2);
attach(pic, legend(pic, black+linewidth(.5pt)),point(pic,NE),15S+15W,Fill(gray(0.98)));
// labels for spin & mass
label(pic, "$M\!=\!"+string(mass)+" M_\odot,\;a\!=\!"+string(spin)+"$", point(pic,NE), 25S+13.6W, fontsize(8pt), NoFill);
label(pic, "$\alpha="+string(alpha)+",\;{\rm inc}\!=\!"+string(inc)+"^\circ$", point(pic,NE), 30S+13W, fontsize(8pt), NoFill);
//label(pic, "$h_f=\rm{BHSPEC}$", point(pic,NE), 35S+14W, fontsize(8pt), NoFill);
label(pic, "$h_f=1.6$", point(pic,NE), 35S+11.5W, fontsize(8pt), NoFill);
*/
return pic;
}
picture P1 = plot_spectrum();
add(P1.fit(),(0,0),N);
import graph;
void read_spectrum(string filename, real[] E, real[] I)
{
file f = input(filename).line();
real[][] raw = f.dimension(0,0);
close(f);
// if row[0] contains only one number, we have module 15 and we skip first 4 lines
// if row[0] contains two numbers, we have module 15i and we skip first 2 lines
int row_offset = (raw[0].length == 1) ? 3 : 2;
real[][] tmp = transpose(raw[row_offset:]);
E.delete();
E.append(tmp[0]*4.135667e-18); // 4.135667e-18 = Hz->keV conversion factor (planck_h/(1e3*electronvolt))
I.delete();
I.append(tmp[1]);
}
// --- settings ----------------
//real mass = 10;
//real spin = 0.65;
//real alpha = 0.100;
//real mdot[] = {0.1, 0.3, 0.6};
//pen colors[] = {black,deepblue,heavyred}; if (settings.gray) colors = new pen[] {black,black,black};
//real inc = 66;
real xmin = 0.01;
real xmax = 30.0;
real ymin = -infinity;
real ymax = infinity;
real size_factor = 1.0;
real size_aspect = 3/2;
defaultpen(currentpen+linewidth(1pt));
usersetting();
pen pen_original = black+linetype("3 3");
pen pen_retrad = black+linetype("0 3");
pen pen_final = black+solid;
//------ plot graph --------------
picture plot_spectrum() {
picture pic;
size(pic, 80mm*size_factor*size_aspect, 80mm*size_factor, IgnoreAspect);
scale(pic, Log, Log);
real[] E;
real[] I;
read_spectrum("../a000/l010/0013/stage1.15", E, I);
draw(pic, graph(pic, E, I), pen_original, "original");
read_spectrum("../a000/l010/0013/stage1.15i", E, I);
draw(pic, graph(pic, E, I), pen_retrad, "retrad");
read_spectrum("../a000/l010/0013/stage2.15", E, I);
draw(pic, graph(pic, E, I), pen_final, "final");
real ymax = 10^(round(log10(max(I)))+0.5);
real ymin = ymax*1e-5;
// label flux curves
// label(pic, rotate(-62)*"$L=0.1$", (log10(9.5),log10(5*ymin)), W, Helvetica()+fontsize(9pt), NoFill);
// label(pic, rotate(-62)*"$L=0.3$", (log10(13.5),log10(5*ymin)), W, Helvetica()+fontsize(9pt), NoFill);
// label(pic, rotate(-62)*"$L=0.6$", (log10(21),log10(5*ymin)), W, Helvetica()+fontsize(9pt), NoFill);
xlimits(pic, xmin, xmax);
ylimits(pic, ymin, ymax);
crop(pic);
xaxis(pic, "energy [keV]", BottomTop, LeftTicks(DefaultFormat), above=true);
yaxis(pic, "flux [erg/cm$^2$/keV/s]", LeftRight, RightTicks, above=true);
/*
// attach legend
legendlinelength = 40;
legendmargin = 6;
Legend l1 = Legend(" NT", black+linewidth(1pt)+dashed+Helvetica()+fontsize(9pt));
Legend l2 = Legend(" SD", black+linewidth(1pt)+solid+Helvetica()+fontsize(9pt));
pic.legend.push(l1);
pic.legend.push(l2);
attach(pic, legend(pic, black+linewidth(.5pt)),point(pic,NE),15S+15W,Fill(gray(0.98)));
// labels for spin & mass
label(pic, "$M\!=\!"+string(mass)+" M_\odot,\;a\!=\!"+string(spin)+"$", point(pic,NE), 25S+13.6W, fontsize(8pt), NoFill);
label(pic, "$\alpha="+string(alpha)+",\;{\rm inc}\!=\!"+string(inc)+"^\circ$", point(pic,NE), 30S+13W, fontsize(8pt), NoFill);
//label(pic, "$h_f=\rm{BHSPEC}$", point(pic,NE), 35S+14W, fontsize(8pt), NoFill);
label(pic, "$h_f=1.6$", point(pic,NE), 35S+11.5W, fontsize(8pt), NoFill);
*/
return pic;
}
picture P1 = plot_spectrum();
add(P1.fit(),(0,0),N);
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment