from IndiceVerre import *
from Dioptres import *
from pylab import *
cat = CatalogueVerre()
n1 = Vide()
n2 = cat.verre['N-BK7']
e = 0.01
c = 1.0
r = 0.05
dioptre1 = Spherique(-e/2,c,n1,n2,r)
dioptre2 = Spherique(e/2,-c,n2,n1,r)
lentille = SystemeCentre()
lentille.ajouter(dioptre1)
lentille.ajouter(dioptre2)
figure()
[x1,z1] = dioptre1.profil(20)
plot(z1,x1,color='b')
[x2,z2] = dioptre2.profil(20)
plot(z2,x2,color='b')
xlabel('z')
ylabel('x')
axis([-0.1,0.1,-0.05,0.05])
grid(True)
            

raie = 'd'
lentille.matriceTransfert(raie)
[Fo,Fi,Ho,Hi] = lentille.foyers(raie)
            

f = FaisceauParalleleMeridien(10,0.0,-1.0,r)
lentille.refractionFaisceau(f,raie)
figure()
plot(z1,x1,color='b')
plot(z2,x2,color='b')
for Rc in f.listeRayons:
    XYZ = Rc.getXYZ(1)
    plot(XYZ[2],XYZ[0],color='y')
xlabel('z') 
ylabel('x')
axis([-1,2,-0.05,0.05])
grid(True)
            

axis([0.96,0.98,-0.001,0.001])
            

def ASL(systeme,raie,r):
    systeme.matriceTransfert(raie)
    [Fo,Fi,Ho,Hi] = systeme.foyers(raie)
    R = Rayon([r,0,-100.0],[0,0,1.0])
    Rc = RayonComplet()
    Rc.ajouter(R)
    systeme.refraction(Rc,raie)
    R2 = Rc.getRayon()
    t = -R2.P[0]/R2.U[0]
    z = R2.P[2]+R2.U[2]*t
    return Fi.z - z
            

asl = ASL(lentille,raie,r)
            

def AslLentille(raie,r,nv,f,c2):
    n1 = Vide()
    c1 = c2+1.0/(f*(nv.indice(raie)-1))
    e = 0.01
    dioptre1 = Spherique(-e/2,c1,n1,nv,r)
    dioptre2 = Spherique(e/2,c2,nv,n1,r)
    lentille = SystemeCentre()
    lentille.ajouter(dioptre1)
    lentille.ajouter(dioptre2)
    return ASL(lentille,raie,r)
            

c2List = []
asl = []
N = 30
for k in range(N):
    c2 = -1.0*k/N
    c2List.append(c2)
    asl.append(AslLentille(raie,r,n2,1.0,c2))
figure()
plot(c2List,asl,color='b')
xlabel('c2')
ylabel('ASL')
grid(True)
axis([-1.0,0.0,0,0.005])
            

c2 = -0.3
f = 1.0
c1 = c2+1.0/(f*(n2.indice(raie)-1))
dioptre1.c = c1
dioptre2.c = c2
f = FaisceauParalleleMeridien(10,0.0,-1.0,r)
lentille.refractionFaisceau(f,raie)
figure()
for Rc in f.listeRayons:
    XYZ = Rc.getXYZ(1)
    plot(XYZ[2],XYZ[0],color='y')
xlabel('z') 
ylabel('x')
axis([0.98,1.01,-0.001,0.001])
grid(True)
            

asl = AslLentille(raie,r,n2,1.0,c2)
            

def ACL(systeme,raie1,raie2):
    systeme.matriceTransfert(raie1)
    [Fo,Fi,Ho,Hi] = systeme.foyers(raie)
    z1 = Fi.z
    systeme.matriceTransfert(raie2)
    [Fo,Fi,Ho,Hi] = systeme.foyers(raie)
    acl = Fi.z - z1
    systeme.matriceTransfert('d')
    [Fo,Fi,Ho,Hi] = systeme.foyers('d')
    fD = Fi.z - Hi.z
    return [acl,fD]
            
 
[acl,fD] = ACL(lentille,'F','C')
            

def AclDoublet(c1,c2,n2,n3,e1,e2,r):
    n1 = Vide()
    dioptre1 = Spherique(0,c1,n1,n2,r)
    dioptre2 = Spherique(e1,c2,n2,n3,r)
    dioptre3 = Spherique(e1+e2,0,n3,n1,r)
    systeme = SystemeCentre()
    systeme.ajouter(dioptre1)
    systeme.ajouter(dioptre2)
    systeme.ajouter(dioptre3)
    [acl,fD] = ACL(systeme,'F','C')
    return [acl,fD]
            

n2 = cat.verre['N-BK7']
n3 = cat.verre['F5']
            

e1 = 0.01
e2 = 0.007
r = 0.05
c1 = 2.4
y = []
c2List = []
N = 10
for k in range(N):
    c2 = -2.0-1.0/N*k
    c2List.append(c2)
    [acl,fD] = AclDoublet(c1,c2,n2,n3,e1,e2,r)
    y.append(acl)
figure()
plot(c2List,y,color='b')
xlabel('c2')
ylabel('ACL')
grid(True)
axis([-2.6,-2.4,-0.001,0.001])
            

c2 = -2.50
dioptre1 = Spherique(0,c1,n1,n2,r)
dioptre2 = Spherique(e1,c2,n2,n3,r)
dioptre3 = Spherique(e1+e2,0,n3,n1,r)
doublet = SystemeCentre()
doublet.ajouter(dioptre1)
doublet.ajouter(dioptre2)
doublet.ajouter(dioptre3)
figure()
[x1,z1] = dioptre1.profil(20)
plot(z1,x1,color='b')
[x2,z2] = dioptre2.profil(20)
plot(z2,x2,color='b')
[x3,z3] = dioptre3.profil(20)
plot(z3,x3,color='b')
xlabel('z')
ylabel('x')
axis([-0.1,0.1,-0.05,0.05])
grid(True)

            

[acl,fD] = ACL(doublet,'F','C')
asl = ASL(doublet,'d',r)
            

[acl,fD] = ACL(doublet,'F','d')
            

figure()
f = FaisceauParalleleMeridien(10,0.0,-1.0,r)
doublet.refractionFaisceau(f,'C')
for Rc in f.listeRayons:
    XYZ = Rc.getXYZ(2)
    plot(XYZ[2],XYZ[0],color='r')
f = FaisceauParalleleMeridien(10,0.0,-1.0,r)
doublet.refractionFaisceau(f,'F')
for Rc in f.listeRayons:
    XYZ = Rc.getXYZ(2)
    plot(XYZ[2],XYZ[0],color='b')
f = FaisceauParalleleMeridien(10,0.0,-1.0,r)
doublet.refractionFaisceau(f,'d')
for Rc in f.listeRayons:
    XYZ = Rc.getXYZ(2)
    plot(XYZ[2],XYZ[0],color='y')
xlabel('z') 
ylabel('x')
axis([0.978,0.982,-0.0001,0.0001])
grid(True)

            

alpha = 0.05
f = FaisceauParalleleCylindre(30,alpha,-0.1,0.02)
doublet.refractionFaisceau(f,'d')
figure()
for Rc in f.listeRayons:
    XYZ = Rc.getXYZ(1.5)
    plot(XYZ[2],XYZ[0],color='y')
xlabel('z') 
ylabel('x')
plot(z1,x1,color='b')
plot(z2,x2,color='b')
plot(z3,x3,color='b')
grid(True)
            

zi = 0.9795
planImage = PlanImage(zi,0.1)
doublet.ajouter(planImage)
def pointsImage(systeme,alpha,r,raie,color):
    f = FaisceauParalleleCylindre(30,alpha,-0.1,r)
    systeme.refractionFaisceau(f,raie)
    x = []
    y = []
    for Rc in f.listeRayons: 
        R = Rc.getRayon()
        x.append(R.P[0])
        y.append(R.P[1])
    plot(x,y,color=color)  
figure(figsize=(7,7))
alpha = 0.005
rList = [0.001,0.005,0.01,0.015,0.02,0.025,0.03,0.035]
for r1 in rList:
    pointsImage(doublet,alpha,r1,'d','y')
xlabel('x')
ylabel('y')
grid(True)
axis([0.004860,0.004880,-0.00001,0.00001])
            

figure(figsize=(7,7))
alpha = 0
for r1 in rList:
    pointsImage(doublet,alpha,r1,'d','y')
xlabel('x')
ylabel('y')
grid(True)
axis([-0.00001,0.00001,-0.00001,0.00001])
            

alpha = 0.02
f = FaisceauParallele(20,alpha,-0.1,0.001)
doublet.refractionFaisceau(f,'d')
figure()
for Rc in f.listeRayons:
    XYZ = Rc.getXYZ(1)
    plot(XYZ[2],XYZ[0],color='y')
xlabel('z')
ylabel('x')
grid(True)
axis([0.978,0.98,0.01948,0.0195])
            

figure()
for Rc in f.listeRayons:
    XYZ = Rc.getXYZ(1)
    plot(XYZ[2],XYZ[1],color='y')
xlabel('z')
ylabel('x')
grid(True)
axis([0.978,0.98,-0.00001,0.00001])
            

planImage.z = 0.9795
alpha = 0.01
f = FaisceauParallele(20,alpha,-0.1,0.05)
doublet.refractionFaisceau(f,'d')
figure(figsize=(8,8))
for Rc in f.listeRayons:
    P = Rc.getRayon().P
    plot([P[0]],[P[1]],color='y',marker='.')
rAiry = 0.000014
angles = linspace(0,2*pi,100)
xi = 0.009748
plot(xi+rAiry*cos(angles),rAiry*sin(angles),color='k')
xlabel('x')
ylabel('y')
grid(True)
axis([0.00969,0.00969+0.0001,-0.00005,0.00005])