9.7.10 Beam tracing sample

C++ code

mglData r, xx, yy, a, im(128), re(128);
const char *ham = "p^2+q^2-x-1+i*0.5*(y+x)*(y>-x)";
r = mglRay(ham, mglPoint(-0.7, -1), mglPoint(0, 0.5), 0.02, 2);
gr->Plot(r.SubData(0), r.SubData(1), "k");
gr->Axis();	gr->Label('x', "\\i x");	gr->Label('y', "\\i z");
// now start beam tracing
re.Fill("exp(-48*x^2)", gr->Min, gr->Max);
a = mglQO2d(ham, re, im, r, 1, 30, &xx, &yy);
gr->CAxis(0, 1);
gr->Dens(xx, yy, a, "wyrRk");
gr->Plot("-x", "k|");
gr->Puts(mglPoint(0, 0.85), "absorption: (x+y)/2 for x+y>0");
gr->Puts(mglPoint(0.7, -0.05), "central ray");
gr->Title("Beam and ray tracing", "C", -1.5);

MGL code

define $1 'p^2+q^2-x-1+i*0.5*(y+x)*(y>-x)'
ray r $1 -0.7 -1 0 0 0.5 0 0.02 2
plot r(0) r(1) 'k'
axis
xlabel '\i x'
ylabel '\i z'
new re 128
new im 128
new xx
new yy
fill re 'exp(-48*x^2)'
qo2d a $1 re im r 1 30 xx yy 
crange 0 1
dens xx yy a 'wyrRk'
fplot '-x' 'k|'
text -0.8 0.85 'absorption: (x+y)/2 for x+y>0' 'L' -1
text 0.5 -0.05 'central ray' 'L' -1
title 'Beam and ray tracing' 'C' -1.5

Pure C code

const char *ham = "p^2+q^2-x-1+i*0.5*(y+x)*(y>-x)";
HMDT r, xx, yy, a, im, re;
im = mgl_create_data_size(128,1,1);
re = mgl_create_data_size(128,1,1);
r = mgl_ray_trace(ham, -0.7, -1., 0., 0., 0.5, 0., 0.02, 2.);
xx = mgl_data_subdata(r,0,-1,-1);
yy = mgl_data_subdata(r,1,-1,-1);
mgl_plot_xy(gr, xx, yy, "k");
mgl_axis(gr,"xyz");
mgl_label(gr,'x',"\\i x");
mgl_label(gr,'y',"\\i z");
mgl_data_fill_eq(gr,re,"exp(-48*x^2)", 0, 0);
a = mgl_qo2d_solve(ham, re, im, r, 1, 30, xx, yy);
mgl_set_caxis(gr,0.,1.);
mgl_dens_xy(gr,xx, yy, a, "wyrRk", -1.);
mgl_fplot(gr,"-x", "k|",100);
mgl_puts(gr, 0., 0.85, 0., "absorption: (x+y)/2 for x+y>0");
mgl_puts(gr, 0.7, -0.05, 0., "central ray");
mgl_title(gr, "Beam and ray tracing", "C", -1.5);
mgl_delete_data(a);     mgl_delete_data(r);
mgl_delete_data(xx);    mgl_delete_data(yy);
mgl_delete_data(im);    mgl_delete_data(re);

Fortran code

integer r, xx, yy, a, im, re, mgl_create_data_size
integer mgl_data_subdata, mgl_ray_trace, mgl_qo2d_solve
character*64 ham
ham = 'p^2+q^2-x-1+i*0.5*(y+x)*(y>-x)'
im = mgl_create_data_size(128,1,1)
re = mgl_create_data_size(128,1,1)
r = mgl_ray_trace(ham, -0.7, -1., 0., 0., 0.5, 0., 0.02, 2.)
xx = mgl_data_subdata(r,0,-1,-1)
yy = mgl_data_subdata(r,1,-1,-1)
call mgl_plot_xy(gr, xx, yy, 'k')
call mgl_axis(gr,'xyz')
call mgl_label(gr,'x','\i x')
call mgl_label(gr,'y','\i z')
call mgl_data_fill_eq(gr,re,'exp(-48*x^2)', 0, 0)
a = mgl_qo2d_solve(ham, re, im, r, 1., 30., xx, yy)
call mgl_set_caxis(gr,0.,1.)
call mgl_dens_xy(gr,xx, yy, a, 'wyrRk', -1.)
call mgl_fplot(gr,'-x', 'k|',100)
call mgl_puts(gr, 0., 0.85, 0., 'absorption: (x+y)/2 for x+y>0')
call mgl_puts(gr, 0.7, -0.05, 0., 'central ray')
call mgl_title(gr, 'Beam and ray tracing', 'C', -1.5)
call mgl_delete_data(a)
call mgl_delete_data(r)
call mgl_delete_data(xx)
call mgl_delete_data(yy)
call mgl_delete_data(im)
call mgl_delete_data(re)

Python

ham = "p^2+q^2-x-1+i*0.5*(y+x)*(y>-x)";
r, xx, yy, a = mglData(), mglData(), mglData(), mglData();
im, re = mglData(128), mglData(128);
r = mglRay(ham, mglPoint(-0.7, -1), mglPoint(0, 0.5), 0.02, 2);
gr.Plot(r.SubData(0), r.SubData(1), "k");
gr.Axis();	gr.Label('x', "\\i x");	gr.Label('y', "\\i z");
gr.Fill(re,"exp(-48*x^2)");
a = mglQO2d(ham, re, im, r, 1, 30, xx, yy);
gr.SetCRange(0, 1);
gr.Dens(xx, yy, a, "wyrRk");
gr.Plot("-x", "k|");
gr.Puts(0, 0.85, 0, "absorption: (x+y)/2 for x+y>0");
gr.Puts(0.7, -0.05, 0, "central ray");
gr.Title("Beam and ray tracing", "C", -1.5);