from math import *
from numpy import *
from pylab import *
xa = []
H = []
for k in range(100):
    x = 1.0/100*k
    xa.append(x)
    H.append(2/pi*(acos(x)-x*sqrt(1-x*x)))
clf()
plot(xa,H)
xlabel("fx/fc")
ylabel("H")
axis([0,2,0,1])
            
from math import *
import Image
from TfdImage import *
from numpy.fft import fft2, ifft2
Nx=400
Ny=400
                

def pupille(r):
    P = zeros((Ny,Nx),dtype=complex)
    nx=Nx/2
    ny=Ny/2
    for k in range(Nx):
        for m in range(Ny):
            if pow((k-nx)*(k-nx)+(m-ny)*(m-ny),0.5) < r:
                P[m,k] = 1.0
    return P
                

D=0.01
np=80
Lx=D*Nx/np
Ly=D*Ny/np
                

P=pupille(np/2)
TFP=matriceCentre(fft2(P))
h=matricePuissance(TFP)
                

Lambda=500e-9
imh = matriceImage(h,3.0,[1.0,1.0,1.0])
xmax = Nx/Lx*0.5*Lambda*1000
figure()
imshow(imh,extent=[-xmax,xmax,-xmax,xmax])
xlabel('X (mrad)')
ylabel('Y (mrad)')
                

H=matriceCentre(fft2(h))
imH=matriceImage(matriceModule(H),1.0,[1.0,1.0,1.0])
figure()
fmax = Nx/(2*xmax)*0.5;
imshow(imH,extent=[-fmax,fmax,-fmax,fmax])
xlabel('fx (1/mrad)')
ylabel('fy (1/mrad)')
                

FTM=imH[Ny/2]
figure()
fx=arange(-fmax,fmax,2*fmax/Nx)
plot(fx,FTM)
xlabel('fx (cycles/mrad)')
ylabel('FTM')
axis([-fmax,fmax,0,1.0])
                

fc = D/Lambda/1000
               

def imageLentille(objet,H):
    TFO = matriceCentre(fft2(objet))
    TFI = zeros((Ny,Nx),dtype=complex)
    for n in range(Nx):
        for l in range(Ny):
            TFI[l,n] = TFO[l,n]*H[l,n]
    Im = matriceCentre(ifft2(matriceBords(TFI)))
    norm=Im.max()
    for n in range(Nx):
        for l in range(Ny):
            Im[l,n]=Im[l,n]/norm
    return Im
                

A=zeros((Ny,Nx),dtype=float)
L0=40.0
ny=Ny/2
for k in range(Nx):
    L=L0-L0/(float(Nx)*1.5)*float(k)
    I=0.5*(1+cos(2*pi*k/L))
    for m in range(ny):
        A[m,k] = I
        if I>0.5:
            A[m+ny,k] = 1.0
figure()
scale = xmax
imshow(matriceImage(A,1.0,[1.0,1.0,1.0]),extent=[-scale,scale,-scale,scale])
xlabel('X (mrad)')
ylabel('Y (mrad)')
                

IA=imageLentille(A,H)
IA=matriceModule(IA)
imIA=matriceImage(IA,1.0,[1.0,1.0,1.0])
figure()
figure()
imshow(imIA,extent=[-scale,scale,-scale,scale])
xlabel('X (mrad)')
ylabel('Y (mrad)')
                

figure()
xlabel('X (mrad)')
ylabel('I')
x=arange(-scale,scale,2*scale/Nx)
plot(x,IA[50])
axis([-scale,scale,0,1])