from CalculImage import * 
cat = CatalogueVerre() 
n1 = Vide() 
n2 = cat.verre['N-BK7'] 
rc = 1033.0
r = 30.0
dioptre1 = Spherique(0,1.0/rc,n1,n2,r)
dioptre2 = Spherique(2,-1.0/rc,n2,n1,r)
sys = SystemeCentre()
sys.ajouterListe([dioptre1,dioptre2])
PObjet = Position(0)
PObjet.infini = True
sys.matriceTransfert('d')
PImage = sys.image(PObjet,'d')
zi = PImage.z
            

zs = dioptre2.z+0.1
R = zi-zs 
sphereImage = SphereImage(zs,[0,0,zi])
sys.ajouter(sphereImage)
planImage = PlanImage(zi,200)
sys.ajouter(planImage)
            

cal = CalculImage(sys)
Nx = Ny = 400
np = 200
rm = 200.0
[P,M,Phi] = cal.surfacePhase(PObjet,0,'d',rm,1,np,Nx)
from TfdImage import *
from numpy.fft import fft2, ifft2
TFP=matriceCentre(fft2(P))
I=matricePuissance(TFP)
Lambda = cal.Lambda['d']
imI = matriceImage(I,2.0,[1.0,1.0,1.0])
figure()
xScale = Nx/(2*rm)*0.5*Lambda*1e-6*R
imshow(imI,extent=[-xScale,xScale,-xScale,xScale])
xlabel('x (mm)')
ylabel('y (mm)') 
            

H = matriceCentre(fft2(I))
imH=matriceImage(matriceModule(H),1.0,[1.0,1.0,1.0])
fxScale = Nx/(2*xScale)*0.5;
imshow(imH,extent=[-fxScale,fxScale,-fxScale,fxScale])
xlabel('fx (cycles/mmm)')
ylabel('fy (cycles/mm)')
            

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

angle = 60
figure(figsize=(6,6))
xlabel('x')
ylabel('y')
xCentre = cal.diagrammeSpots(PObjet,angle,['d'],0,20)
rAiry = cal.airy(PObjet,'d')
grid(True)
scale=0.1
axis([-scale,scale,-scale,scale])
            

sphereImage.PC = [xCentre,0,zi]
Nx = Ny = 400
np = 200
rm = 200.0
[P,M,Phi] = cal.surfacePhase(PObjet,angle,'d',rm,1,np,Nx)
clf()
xlabel('x')
ylabel('y')
title('Phase')
imshow(Phi)
colorbar()

            

TFP=matriceCentre(fft2(P))
I=matricePuissance(TFP)
Lambda = cal.Lambda['d']
imI = matriceImage(I,2.0,[1.0,1.0,1.0])
figure()
xScale = Nx/(2*rm)*0.5*Lambda*1e-6*R
imshow(imI,extent=[-xScale,xScale,-xScale,xScale])
xlabel('x (mm)')
ylabel('y (mm)') 
            

H = matriceCentre(fft2(I))
imH=matriceImage(matriceModule(H),1.0,[1.0,1.0,1.0])
fxScale = Nx/(2*xScale)*0.5;
imshow(imH,extent=[-fxScale,fxScale,-fxScale,fxScale])
xlabel('fx (cycles/mmm)')
ylabel('fy (cycles/mm)')
            

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