PSF eigenvalue

The objective is to calculate eigenvalue of a PSF matrix.


In [1]:
using Revise


┌ Info: Precompiling Revise [295af30f-e4ad-537b-8983-00126c2a3abe]
└ @ Base loading.jl:1273

In [2]:
using PyPlot, LinearAlgebra
pygui(false)


┌ Info: Precompiling PyPlot [d330b81b-6aea-500a-939a-2ce795aea3ee]
└ @ Base loading.jl:1273
Out[2]:
false

In [3]:
include("psfeig.jl")
include("utils.jl")


Out[3]:
Main.utils

Reading Modeled PSFs


In [96]:
path = "/home/bruno/Work/Luana_PSF/psfs_su5_size71/all_psfs_su5_size71_com_ph_interp.rsf@"
image = utils.read_bin(path,960,5300);

In [97]:
utils.sisshow(image, vmin=-0.1, vmax=0.1)



In [98]:
b1=70; b2=70; k1=36; k2=36; l1=920; l2=5260;
println("psf grid: ($(length(k1:b1:l1)),$(length(k2:b2:l2)))")
psfs = psfeig.psf_chop(image,b1,b2,k1,k2,l1,l2);


psf grid: (13,75)

In [99]:
size(psfs.val,3)


Out[99]:
975

In [100]:
utils.sisshow(psfs.val[:,:,9*13+1], psfs.val[:,:,9*13+5], psfs.val[:,:,9*13+9], psfs.val[:,:,9*13+13], perc=99.9)



In [11]:
o=36; l=12
imshow(psfs.val[(o-l):(o+l),(o-l):(o+l),1], cmap="Greys")
grid()



In [12]:
size(psfs.val[(o-l):(o+l),(o-l):(o+l),1])


Out[12]:
(25, 25)

In [13]:
findmax(psfs.val[(o-l):(o+l),(o-l):(o+l),1])


Out[13]:
(3.3668396f0, CartesianIndex(13, 13))

Eigenvalues of Modeled PSFs


In [101]:
function plot_psf_with_condition_number(ipsf)
    P = psfs.val[(o-l):(o+l),(o-l):(o+l),ipsf];
    G = psfeig.psf_to_matrix(n,P);
    figure()
    imshow(psfs.val[(o-l):(o+l),(o-l):(o+l),ipsf],
            interpolation="bicubic", cmap="Greys", vmin=-1, vmax=1)
    colorbar()
    grid()
    title("PSF #$(ipsf): condition number = $(cond(G))")  
end


Out[101]:
plot_psf_with_condition_number (generic function with 1 method)

Distict PSFs


In [15]:
o=36; l=12; n=25; 

for ipsf in [9*13+1, 9*13+5, 9*13+9, 9*13+13]
    plot_psf_with_condition_number(ipsf)
end



In [17]:
o=36; l=12; n=25;

for ipsf in [9*13+1, 10*13+1, 9*13+2, 10*13+2]
    plot_psf_with_condition_number(ipsf)
end


Eigenvalues of 4x4 PSFs

Distinct PSFs


In [102]:
o=36; l=12; n = 100

P1 = psfs.val[(o-l):(o+l),(o-l):(o+l), 9*13+1];
P2 = psfs.val[(o-l):(o+l),(o-l):(o+l), 9*13+5];
P3 = psfs.val[(o-l):(o+l),(o-l):(o+l), 9*13+9];
P4 = psfs.val[(o-l):(o+l),(o-l):(o+l), 9*13+13];

In [103]:
G = psfeig.psf_to_matrix(n,P1,P2,P3,P4); varinfo(r"G")


Out[103]:
name size summary
G 762.939 MiB 10000×10000 Array{Float64,2}

In [27]:
cond(G)


Out[27]:
3.152389269824061e8

In [28]:
imshow(reshape(G[(75-1)*100+75,:],100,100)); grid()


Neighbour PSFs


In [56]:
o=36; l=12; n = 100

P1 = psfs.val[(o-l):(o+l),(o-l):(o+l),  9*13+1];
P2 = psfs.val[(o-l):(o+l),(o-l):(o+l), 10*13+1];
P3 = psfs.val[(o-l):(o+l),(o-l):(o+l),  9*13+2];
P4 = psfs.val[(o-l):(o+l),(o-l):(o+l), 10*13+2];

In [57]:
G = psfeig.psf_to_matrix(n,P1,P2,P3,P4); varinfo(r"G")


Out[57]:
name size summary
G 762.939 MiB 10000×10000 Array{Float64,2}

In [31]:
cond(G)


Out[31]:
1.8344483898310825e7

In [32]:
imshow(reshape(G[(50-1)*100+50,:],100,100)); grid()


PSF Eigenvalue Developmment


In [39]:
n = 25;
P = psfeig.bluker(12,1.0);

In [107]:
G = psfeig.psf_to_matrix(n,P); varinfo(r"G")
G2 = psfeig.psf_to_matrix(n,P,P,P,P); varinfo(r"G2")
G3 = psfeig.psf_to_matrix(n,P,0.5*P,0.2*P,0.1*P); varinfo(r"G3")


Out[107]:
name size summary
G3 2.980 MiB 625×625 Array{Float64,2}

In [109]:
cond(G)


Out[109]:
4167.911018197733

In [110]:
cond(G2)


Out[110]:
4167.911018197979

In [111]:
cond(G3)


Out[111]:
11380.202021198646

In [105]:
imshow(reshape(G[(11-1)*25+11,:],25,25)); grid()



In [103]:
imshow(reshape(G4[(11-1)*25+11,:],25,25)); grid()



In [40]:
imshow(P)


Out[40]:
PyObject <matplotlib.image.AxesImage object at 0x7fbe251ea410>

In [41]:
findmax(P)


Out[41]:
(0.15915494309189535, CartesianIndex(13, 13))

In [42]:
size(P)


Out[42]:
(25, 25)

Eigenvalue Timing

$n$ Time Size
64 27s 128Mb
96 162s 324Mb

In [36]:
n=100; x = rand(Float32, n*n,n*n); varinfo(r"x")


Out[36]:
name size summary
x 381.470 MiB 10000×10000 Array{Float32,2}

In [83]:
@time v = eigvals(x);


165.460979 seconds (899.31 k allocations: 428.171 MiB, 0.01% gc time)

In [84]:
x = nothing

In [85]:
GC.gc()

Evaluate the condition number of blur matrix


In [35]:
N = 50;
k = 12;

nn = 11;
c = zeros(nn);
sigma = range(0.1, stop=2.0, length=nn);
@time for i = 1:nn
    P = psfeig.bluker(k,sigma[i]);
    G = psfeig.psf_to_matrix(N,P);
    c[i] = cond(G);
end

semilogy(sigma,c);


 47.223515 seconds (81.13 k allocations: 1.044 GiB, 0.16% gc time)

Gulfaks Reflectivity


In [22]:
using FileIO

In [53]:
gulfacks = load(File(format"PNG", "gulfacks.png"));
refl = [ Float32(gulfacks[i,j].r) for i=1:100,j=1:100];

In [43]:
imshow([refl[i,j].r for i=1:100,j=1:100], cmap="Greys")
colorbar()


Out[43]:
PyObject <matplotlib.colorbar.Colorbar object at 0x7f3a8a433f10>

In [62]:
function hprod(G,r)
    n1,n2 = size(r)
    m1,m2 = size(G)
    @assert(m1 == n1*n2, "Hessian and reflectivity are not compatible")
    return reshape(G*reshape(r,n1*n2),n1,n2)
end

function hinv(G,r)
    n1,n2 = size(r)
    m1,m2 = size(G)
    @assert(m1 == n1*n2, "Hessian and reflectivity are not compatible")
    return reshape(G \ reshape(r,n1*n2),n1,n2)
end


Out[62]:
hinv (generic function with 1 method)

In [69]:
hr = hprod(G,refl);

In [55]:
imshow(hr, cmap="Greys")


Out[55]:
PyObject <matplotlib.image.AxesImage object at 0x7f3a8a6fbfd0>

In [105]:
hr = hprod(G,refl);

In [106]:
imshow(hr, cmap="Greys"); title("Com PH")


Out[106]:
PyObject Text(0.5, 1, 'Com PH')

In [70]:
refl_inv = hinv(G,hr);
imshow(refl_inv, cmap="Greys"); colorbar()


Out[70]:
PyObject <matplotlib.colorbar.Colorbar object at 0x7f3a8a304510>

In [58]:
hr = hprod(G,refl);
imshow(hr, cmap="Greys")


Out[58]:
PyObject <matplotlib.image.AxesImage object at 0x7f3a8a4422d0>

In [63]:
refl_inv = hinv(G,hr);
imshow(refl_inv, cmap="Greys"); colorbar();


Out[63]:
100×100 Array{Float64,2}:
 0.12549    1.48024e-11  -3.15238e-11  …   1.73618e-12  -1.27932e-12  0.12549
 0.12549   -3.07931e-11   6.98919e-11     -1.797e-12     9.49101e-13  0.12549
 0.12549    2.7429e-11   -6.22348e-11     -1.13687e-13   2.31988e-12  0.12549
 0.12549   -6.03352e-12  -3.23145e-12     -1.16895e-12  -4.61284e-12  0.12549
 0.12549   -4.24967e-12   7.16542e-11      3.14014e-12   5.59086e-12  0.12549
 0.12549   -1.79864e-11  -8.82675e-11  …  -7.41961e-12  -2.22851e-12  0.12549
 0.12549    8.68955e-11   9.75597e-12      4.20678e-12   4.44889e-14  0.12549
 0.12549   -1.57936e-10   9.27099e-11      1.15617e-12   1.33211e-12  0.12549
 0.12549    1.69316e-10  -1.03901e-10     -3.32198e-12  -1.86793e-12  0.12549
 0.12549   -1.18314e-10   1.58127e-11      2.36034e-13   2.06915e-12  0.12549
 0.12549    5.25272e-11   8.37207e-11  …   4.20772e-12  -1.59966e-12  0.12549
 0.12549   -1.71067e-11  -1.14184e-10     -6.94803e-12   1.66057e-12  0.12549
 0.12549    2.64775e-11   5.17066e-11      2.75835e-12  -1.10247e-12  0.12549
 ⋮                                     ⋱                                     
 0.12549   -3.71283e-11   1.14502e-10      3.60742e-12   9.27299e-14  0.12549
 0.12549    2.17572e-11  -9.68686e-11      4.51247e-12  -2.14178e-12  0.12549
 0.12549    7.33993e-12   6.18708e-11  …  -6.69582e-12   4.22414e-12  0.12549
 0.223529   0.47451       0.513726        -5.58872e-13  -3.67194e-12  0.12549
 0.439216   0.8           0.870588         6.03539e-12   3.92719e-12  0.12549
 0.701961   0.196078      0.121569        -1.23207e-11  -1.49318e-12  0.12549
 0.470588   0.0352941     2.40384e-11      9.88349e-12   1.63719e-12  0.12549
 0.12549    5.77961e-11  -4.95682e-11  …  -9.26668e-12  -4.51412e-14  0.12549
 0.12549   -9.47334e-11   8.3644e-11       9.7804e-12   -1.22203e-12  0.12549
 0.141176   0.2           0.203922        -1.23869e-11   3.18353e-12  0.12549
 0.435294   0.858824      0.929412         0.12549       0.12549      0.25098
 1.0        1.0           1.0              1.0           1.0          1.0    

In [66]:



Out[66]:
PyObject <matplotlib.colorbar.Colorbar object at 0x7f3a8a19f6d0>

In [72]:
using IterativeSolvers


┌ Info: Precompiling IterativeSolvers [42fd0dbc-a981-5370-80f2-aaf504508153]
└ @ Base loading.jl:1273

In [90]:
refl_inv = lsqr(G, hr[:], maxiter=100, verbose=true);


=== lsqr ===
iter	resnorm		  anorm		  cnorm		  rnorm
  1	8.04e+01	9.47e-01	2.69e-01	6.60e-01
  2	6.30e+01	3.96e-01	1.04e-01	5.17e-01
  3	5.13e+01	3.21e-01	6.46e-02	4.21e-01
  4	4.23e+01	2.50e-01	4.55e-02	3.47e-01
  5	3.68e+01	2.05e-01	3.48e-02	3.02e-01
  6	3.18e+01	1.66e-01	2.79e-02	2.61e-01
  7	2.81e+01	1.36e-01	2.35e-02	2.31e-01
  8	2.50e+01	1.20e-01	2.03e-02	2.05e-01
  9	2.20e+01	1.09e-01	1.79e-02	1.81e-01
 10	1.94e+01	9.91e-02	1.60e-02	1.59e-01
 11	1.75e+01	8.76e-02	1.45e-02	1.44e-01
 12	1.62e+01	9.03e-02	1.31e-02	1.33e-01
 13	1.50e+01	8.87e-02	1.19e-02	1.23e-01
 14	1.40e+01	7.57e-02	1.09e-02	1.15e-01
 15	1.30e+01	6.52e-02	1.00e-02	1.07e-01
 16	1.23e+01	6.40e-02	9.34e-03	1.01e-01
 17	1.14e+01	5.74e-02	8.70e-03	9.35e-02
 18	1.07e+01	5.74e-02	8.17e-03	8.79e-02
 19	1.02e+01	5.40e-02	7.68e-03	8.34e-02
 20	9.59e+00	4.53e-02	7.25e-03	7.87e-02
 21	9.11e+00	4.52e-02	6.87e-03	7.48e-02
 22	8.67e+00	4.23e-02	6.53e-03	7.11e-02
 23	8.27e+00	3.85e-02	6.22e-03	6.79e-02
 24	7.92e+00	3.83e-02	5.93e-03	6.50e-02
 25	7.60e+00	3.62e-02	5.67e-03	6.24e-02
 26	7.33e+00	3.59e-02	5.42e-03	6.02e-02
 27	7.08e+00	3.24e-02	5.19e-03	5.81e-02
 28	6.83e+00	3.13e-02	4.98e-03	5.61e-02
 29	6.60e+00	3.04e-02	4.79e-03	5.42e-02
 30	6.38e+00	2.76e-02	4.61e-03	5.24e-02
 31	6.20e+00	2.86e-02	4.44e-03	5.09e-02
 32	6.03e+00	2.63e-02	4.28e-03	4.95e-02
 33	5.88e+00	2.71e-02	4.13e-03	4.83e-02
 34	5.71e+00	2.52e-02	3.98e-03	4.69e-02
 35	5.58e+00	2.74e-02	3.85e-03	4.58e-02
 36	5.45e+00	2.61e-02	3.71e-03	4.48e-02
 37	5.32e+00	2.36e-02	3.59e-03	4.36e-02
 38	5.19e+00	2.27e-02	3.48e-03	4.26e-02
 39	5.07e+00	2.12e-02	3.38e-03	4.16e-02
 40	4.95e+00	2.16e-02	3.28e-03	4.06e-02
 41	4.85e+00	2.16e-02	3.18e-03	3.98e-02
 42	4.75e+00	1.99e-02	3.09e-03	3.90e-02
 43	4.63e+00	2.06e-02	3.00e-03	3.80e-02
 44	4.53e+00	2.23e-02	2.92e-03	3.72e-02
 45	4.47e+00	2.49e-02	2.83e-03	3.67e-02
 46	4.39e+00	2.05e-02	2.74e-03	3.60e-02
 47	4.29e+00	1.89e-02	2.67e-03	3.52e-02
 48	4.20e+00	1.93e-02	2.60e-03	3.45e-02
 49	4.11e+00	1.83e-02	2.53e-03	3.37e-02
 50	4.02e+00	2.12e-02	2.47e-03	3.30e-02
 51	3.97e+00	2.18e-02	2.40e-03	3.26e-02
 52	3.90e+00	1.76e-02	2.34e-03	3.20e-02
 53	3.82e+00	1.66e-02	2.29e-03	3.13e-02
 54	3.73e+00	1.65e-02	2.24e-03	3.06e-02
 55	3.66e+00	1.69e-02	2.19e-03	3.00e-02
 56	3.59e+00	1.69e-02	2.14e-03	2.95e-02
 57	3.52e+00	1.64e-02	2.10e-03	2.89e-02
 58	3.45e+00	1.63e-02	2.05e-03	2.83e-02
 59	3.38e+00	1.56e-02	2.01e-03	2.77e-02
 60	3.30e+00	1.61e-02	1.97e-03	2.71e-02
 61	3.24e+00	1.83e-02	1.93e-03	2.66e-02
 62	3.20e+00	1.68e-02	1.89e-03	2.63e-02
 63	3.15e+00	1.53e-02	1.85e-03	2.58e-02
 64	3.09e+00	1.52e-02	1.82e-03	2.54e-02
 65	3.02e+00	1.51e-02	1.78e-03	2.48e-02
 66	2.97e+00	1.58e-02	1.75e-03	2.44e-02
 67	2.93e+00	1.68e-02	1.72e-03	2.41e-02
 68	2.89e+00	1.45e-02	1.68e-03	2.37e-02
 69	2.84e+00	1.38e-02	1.65e-03	2.33e-02
 70	2.79e+00	1.39e-02	1.63e-03	2.29e-02
 71	2.74e+00	1.30e-02	1.60e-03	2.25e-02
 72	2.69e+00	1.29e-02	1.57e-03	2.21e-02
 73	2.64e+00	1.35e-02	1.55e-03	2.17e-02
 74	2.59e+00	1.33e-02	1.53e-03	2.13e-02
 75	2.55e+00	1.33e-02	1.50e-03	2.09e-02
 76	2.51e+00	1.57e-02	1.48e-03	2.06e-02
 77	2.48e+00	1.43e-02	1.45e-03	2.04e-02
 78	2.45e+00	1.58e-02	1.43e-03	2.01e-02
 79	2.43e+00	1.27e-02	1.40e-03	1.99e-02
 80	2.39e+00	1.27e-02	1.38e-03	1.96e-02
 81	2.35e+00	1.11e-02	1.36e-03	1.93e-02
 82	2.32e+00	1.15e-02	1.34e-03	1.90e-02
 83	2.28e+00	1.25e-02	1.32e-03	1.88e-02
 84	2.27e+00	1.42e-02	1.30e-03	1.86e-02
 85	2.24e+00	1.10e-02	1.28e-03	1.84e-02
 86	2.21e+00	1.27e-02	1.26e-03	1.81e-02
 87	2.19e+00	1.09e-02	1.24e-03	1.80e-02
 88	2.16e+00	1.01e-02	1.23e-03	1.77e-02
 89	2.13e+00	1.01e-02	1.21e-03	1.75e-02
 90	2.10e+00	1.02e-02	1.20e-03	1.72e-02
 91	2.07e+00	1.03e-02	1.18e-03	1.70e-02
 92	2.04e+00	1.05e-02	1.17e-03	1.68e-02
 93	2.02e+00	1.27e-02	1.15e-03	1.66e-02
 94	2.00e+00	1.02e-02	1.14e-03	1.64e-02
 95	1.97e+00	1.02e-02	1.12e-03	1.62e-02
 96	1.94e+00	1.04e-02	1.11e-03	1.60e-02
 97	1.92e+00	1.12e-02	1.09e-03	1.58e-02
 98	1.90e+00	1.05e-02	1.08e-03	1.56e-02
 99	1.88e+00	9.65e-03	1.07e-03	1.54e-02
100	1.86e+00	1.02e-02	1.05e-03	1.52e-02


In [107]:
refl_inv = lsqr(G, hr[:], maxiter=100, verbose=true);


=== lsqr ===
iter	resnorm		  anorm		  cnorm		  rnorm
  1	1.27e+02	5.93e-01	2.19e-01	5.97e-01
  2	8.75e+01	4.30e-01	1.03e-01	4.13e-01
  3	6.14e+01	2.84e-01	6.66e-02	2.90e-01
  4	4.86e+01	2.36e-01	4.91e-02	2.29e-01
  5	3.88e+01	2.07e-01	3.81e-02	1.83e-01
  6	3.30e+01	2.02e-01	3.07e-02	1.56e-01
  7	2.88e+01	1.45e-01	2.52e-02	1.36e-01
  8	2.51e+01	1.27e-01	2.16e-02	1.18e-01
  9	2.23e+01	1.12e-01	1.89e-02	1.05e-01
 10	2.01e+01	1.02e-01	1.67e-02	9.47e-02
 11	1.80e+01	8.78e-02	1.50e-02	8.48e-02
 12	1.66e+01	8.23e-02	1.36e-02	7.82e-02
 13	1.53e+01	7.48e-02	1.24e-02	7.21e-02
 14	1.42e+01	7.01e-02	1.14e-02	6.68e-02
 15	1.34e+01	7.28e-02	1.05e-02	6.31e-02
 16	1.27e+01	6.60e-02	9.63e-03	5.97e-02
 17	1.19e+01	5.80e-02	8.92e-03	5.60e-02
 18	1.11e+01	5.23e-02	8.36e-03	5.25e-02
 19	1.05e+01	5.18e-02	7.85e-03	4.96e-02
 20	9.99e+00	5.08e-02	7.39e-03	4.71e-02
 21	9.53e+00	4.63e-02	6.97e-03	4.50e-02
 22	9.06e+00	4.48e-02	6.59e-03	4.28e-02
 23	8.65e+00	4.18e-02	6.26e-03	4.08e-02
 24	8.23e+00	3.98e-02	5.95e-03	3.88e-02
 25	7.86e+00	3.93e-02	5.68e-03	3.71e-02
 26	7.55e+00	3.68e-02	5.43e-03	3.56e-02
 27	7.27e+00	3.71e-02	5.18e-03	3.43e-02
 28	6.99e+00	3.53e-02	4.96e-03	3.30e-02
 29	6.72e+00	3.55e-02	4.75e-03	3.17e-02
 30	6.46e+00	3.25e-02	4.57e-03	3.05e-02
 31	6.19e+00	3.37e-02	4.39e-03	2.92e-02
 32	5.96e+00	3.26e-02	4.22e-03	2.81e-02
 33	5.73e+00	3.02e-02	4.07e-03	2.70e-02
 34	5.51e+00	3.11e-02	3.93e-03	2.60e-02
 35	5.31e+00	3.03e-02	3.80e-03	2.51e-02
 36	5.11e+00	3.12e-02	3.66e-03	2.41e-02
 37	4.93e+00	3.10e-02	3.54e-03	2.32e-02
 38	4.73e+00	3.08e-02	3.43e-03	2.23e-02
 39	4.55e+00	2.97e-02	3.32e-03	2.15e-02
 40	4.38e+00	3.19e-02	3.22e-03	2.07e-02
 41	4.22e+00	2.83e-02	3.12e-03	1.99e-02
 42	4.04e+00	2.67e-02	3.03e-03	1.91e-02
 43	3.88e+00	2.70e-02	2.95e-03	1.83e-02
 44	3.75e+00	2.69e-02	2.87e-03	1.77e-02
 45	3.62e+00	2.77e-02	2.80e-03	1.71e-02
 46	3.51e+00	2.66e-02	2.73e-03	1.66e-02
 47	3.39e+00	2.43e-02	2.66e-03	1.60e-02
 48	3.30e+00	2.71e-02	2.59e-03	1.56e-02
 49	3.23e+00	2.58e-02	2.52e-03	1.52e-02
 50	3.15e+00	2.44e-02	2.46e-03	1.49e-02
 51	3.06e+00	2.13e-02	2.40e-03	1.45e-02
 52	2.98e+00	2.00e-02	2.35e-03	1.41e-02
 53	2.91e+00	2.02e-02	2.30e-03	1.37e-02
 54	2.84e+00	1.96e-02	2.25e-03	1.34e-02
 55	2.78e+00	1.91e-02	2.21e-03	1.31e-02
 56	2.72e+00	1.84e-02	2.16e-03	1.28e-02
 57	2.67e+00	2.05e-02	2.12e-03	1.26e-02
 58	2.64e+00	1.82e-02	2.07e-03	1.24e-02
 59	2.59e+00	1.69e-02	2.03e-03	1.22e-02
 60	2.54e+00	1.59e-02	1.99e-03	1.20e-02
 61	2.50e+00	1.49e-02	1.96e-03	1.18e-02
 62	2.45e+00	1.51e-02	1.92e-03	1.16e-02
 63	2.42e+00	1.44e-02	1.89e-03	1.14e-02
 64	2.38e+00	1.42e-02	1.86e-03	1.12e-02
 65	2.34e+00	1.36e-02	1.82e-03	1.11e-02
 66	2.31e+00	1.33e-02	1.79e-03	1.09e-02
 67	2.28e+00	1.32e-02	1.76e-03	1.07e-02
 68	2.25e+00	1.27e-02	1.73e-03	1.06e-02
 69	2.22e+00	1.24e-02	1.71e-03	1.05e-02
 70	2.19e+00	1.19e-02	1.68e-03	1.03e-02
 71	2.16e+00	1.27e-02	1.65e-03	1.02e-02
 72	2.14e+00	1.24e-02	1.62e-03	1.01e-02
 73	2.11e+00	1.16e-02	1.60e-03	9.97e-03
 74	2.09e+00	1.16e-02	1.57e-03	9.86e-03
 75	2.06e+00	1.19e-02	1.55e-03	9.74e-03
 76	2.04e+00	1.28e-02	1.52e-03	9.65e-03
 77	2.02e+00	1.19e-02	1.50e-03	9.55e-03
 78	2.00e+00	1.21e-02	1.47e-03	9.44e-03
 79	1.98e+00	1.11e-02	1.45e-03	9.34e-03
 80	1.95e+00	1.14e-02	1.43e-03	9.22e-03
 81	1.93e+00	1.16e-02	1.40e-03	9.11e-03
 82	1.91e+00	1.22e-02	1.38e-03	9.00e-03
 83	1.88e+00	1.28e-02	1.36e-03	8.89e-03
 84	1.86e+00	1.19e-02	1.34e-03	8.80e-03
 85	1.84e+00	1.18e-02	1.32e-03	8.69e-03
 86	1.82e+00	1.23e-02	1.30e-03	8.59e-03
 87	1.80e+00	1.26e-02	1.27e-03	8.48e-03
 88	1.78e+00	1.44e-02	1.25e-03	8.39e-03
 89	1.76e+00	1.12e-02	1.23e-03	8.30e-03
 90	1.74e+00	1.10e-02	1.21e-03	8.19e-03
 91	1.71e+00	1.09e-02	1.20e-03	8.08e-03
 92	1.69e+00	1.13e-02	1.18e-03	7.98e-03
 93	1.67e+00	1.03e-02	1.17e-03	7.88e-03
 94	1.65e+00	1.03e-02	1.15e-03	7.78e-03
 95	1.63e+00	1.07e-02	1.14e-03	7.68e-03
 96	1.61e+00	1.17e-02	1.12e-03	7.60e-03
 97	1.60e+00	1.16e-02	1.10e-03	7.54e-03
 98	1.58e+00	1.01e-02	1.09e-03	7.47e-03
 99	1.56e+00	9.48e-03	1.07e-03	7.37e-03
100	1.55e+00	9.97e-03	1.06e-03	7.29e-03


In [108]:
imshow(reshape(refl_inv,100,100), cmap="Greys", vmin=0, vmax=1); colorbar(); title("Com PH")


Out[108]:
PyObject Text(0.5, 1, 'Com PH')

In [93]:
imshow(reshape(refl_inv,100,100), cmap="Greys", vmin=0, vmax=1); colorbar();



In [ ]: