# -*- coding: utf-8 -*-
# \"""
# \Copyright
# \
# \Rewrite by Alessandro Mirone from the original code by
# \(c) Francesca Mastropietro ESRF-CEA/INAC
# \Wave Front Characterization
# \Francesca Mastropietro/Ana Diaz/Virginie Chamard
# \Simulation of wavefront after FZP and propagation to the detector for the coherent illuminated Fresnl Zone Plate
# \
# \"""
import numpy
import math
import cPickle
from numpy import fft
import h5py
import sys
###### PARAMETERS ######
"""
Energy keV, wavelength m, wavevector m**-1, FZP diameter m, outermost zone m, focal distance m
plane xy ==> FZP plane; z: direction of propagation
"""
import math
# N = 2^a*3^b*5^c*7^d*11^e*13^f, where e + f = 0 or 1, and the rest
# of the exponents are arbitrary
def roundfft(N):
MA,MB,MC,MD,ME,MF = 0,0,0,0,0,0
FA,FB,FC,FD,FE,FFF = 2,3,5,7,11,13
DIFF=9999999999
RES=1
R0=1
AA=1
for A in range( int(math.log(N)/math.log(FA)+2) ):
BB=AA
for B in range( int(math.log(N)/math.log(FB)+2) ):
CC=BB
for C in range( int(math.log(N)/math.log(FC)+2) ):
DD=CC
for D in range( int(math.log(N)/math.log(FD)+2) ):
EE=DD
for E in range( 2 ):
FF=EE
for F in range( 2 -E ):
if FF>=N and DIFF> abs(N-FF):
MA,MB,MC,MD,ME,MF = A,B,C,D,E,F
DIFF=abs(N-FF)
RES=FF
if FF > N: break
FF = FF*FFF
if EE > N: break
EE = EE*FE
if DD > N: break
DD = DD*FD
if CC > N: break
CC = CC*FC
if BB > N: break
BB = BB*FB
if AA > N: break
AA = AA*FA
return RES
[docs]class Parameters :
""" The input variables are read from the input file. The input file is written with python syntax.
The input file has to set the following variables, examples in yellow
* E ::
E=8.0
# units are KeV
* diam ::
diam=200*10**(-6)
# units are meters
* delta_min : the size of the smallest FZP feature ::
delta_min=70*10**(-9)
**All distance units are meters**
* OSA : pinhole diameter ::
OSA=50*10**(-6)
* CS : diameter of central stop ::
CS=65*10**(-6)
* phase : a non essential phase ::
phase = 0.0
* off_set : slit offset horizontally. A typical choice is (CS/2+slit_h/2).
off_set= 50.0*10**(-6)
* slit_h : horizonal slit aperture ::
slit_h=20e-6
* slit_v : vertical slit aperture ::
slit_v=20e-6
* slit_far : slit-FZP distance ::
slit_far= 0.2
* z_fd : focus-detector distance
z_fd= 4.0
* N : the number of sampling point for each dimension. Must be mor or less adapted to the fzp finest feature ::
N=10000
* FILE_FOR_FOCAL_DEPTH_PROFILING , on this file the program will output a 3D profile of the focal region.
if set to None it will not be calculated ::
FILE_FOR_FOCAL_DEPTH_PROFILING = "focal_region"
* range_x_focal_dep and range_z_focal_dep, Nz_focal_dep : in case the above file name is not None, they give the sizes
of the 3D box where the focal intensity is calculated, and the number of slices taken along the longitudinal axis ::
range_x_focal_dep=6.2e-6
range_z_focal_dep=.5*10**(-2)
Nz_focal_dep = 200
* FILE_FOR_DETECTOR : the field at the detector will be written on this file ::
FILE_FOR_DETECTOR ="detector_image"
* detector_pixel_size : the variable name speaks on its own. This determines also the width of the calculated wavefront ::
detector_pixel_size=5.0e-6
"""
if(sys.argv[0][-12:]!="sphinx-build"):
exec(open(sys.argv[1],"r").read())
print "----- READAPTING N for fft "
print "..... OLD:",N,
N=roundfft(N)
print "...... NEW:",N
allowed_pars = ['CS', 'E', 'FILE_FOR_DETECTOR', 'FILE_FOR_FOCAL_DEPTH_PROFILING', 'N', 'Nz_focal_dep', 'OSA', 'delta_min', 'diam', 'off_set', 'phase', 'range_x_focal_dep', 'range_z_focal_dep', 'slit_far', 'slit_h', 'slit_v', 'z_fd', 'detector_pixel_size' ]
if ( (set(locals().keys())-set(allowed_pars+["allowed_pars","__module__"]) )!=set([]) ):
print " PROBLEM: extra parameters have been set by input :", set(dir())-set(allowed_pars+["allowed_pars","__module__"])
raise Exception,""
dic={}
exec(open(sys.argv[1],"r").read() ,globals(),dic)
print " ########################################################################################## "
print " ### Now printing parameters. ################# "
print " "
for key in allowed_pars:
if key in dic.keys():
print key," = " , locals()[key]
else:
print key," = " , "NOT SET !!!!" ," ***********"
print " \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\n"
@staticmethod
[docs] def get_derived_pars():
""" **Derived Parameters**
Besides the parameters set by input, there are other very important
Parameters which are deduced from the input ones
* lambd=12.398/E
* k=2*math.pi/lambd
* focus=(diam*delta_min/lambd
* OSA_position = focus*5.0/6.0
* step=diam/N
* range_x_focus2dect = (12.398/E)*10**(-10)*z_fd/(detector_pixel_size
"""
class pippo:
pass
self=Parameters()
res=pippo()
res.lambd=12.398/self.E*10**(-10)
res.k=2*math.pi/res.lambd
res.focus=(self.diam*self.delta_min/res.lambd)
res.OSA_position = res.focus*5.0/6.0
print " focus " , res.focus
res.step=self.diam/self.N
res.range_x_focus2dect = (12.398/self.E)*10**(-10)*self.z_fd/(self.detector_pixel_size)
return res
#############################################################
# construction of the Fresnel lens
# ----------------------------------
def build_FZ(Parameters, derived_pars , xs ):
radius2= (xs*xs)[:,None] +(xs*xs)[None,:]
path_function = radius2 /(2.0*derived_pars.focus)
path_function = numpy.mod(path_function, derived_pars.lambd)
FZP_2D = numpy.less( derived_pars.lambd/2.0, path_function )
maxradius2 = (Parameters.diam/2.0)**2
FZP_2D = FZP_2D*( numpy.less(radius2,maxradius2 ) ) + numpy.less(maxradius2,radius2 )
CS_2D = radius2 > (Parameters.CS /2.0 )**2
OSA_2D = radius2 < (Parameters.OSA/2.0 )**2
return { "FZP_2D": FZP_2D ,"CS_2D": CS_2D ,"OSA_2D": OSA_2D }
## ----------------------------------------------------------------
#######################################################################
def scrivi_linea(X ,M, nome , nlinea, dire) :
if dire==0:
Y = M[:,nlinea]
elif dire==1:
Y = M[nlinea,:]
else:
raise Exception
f=open(nome,"w")
for x,y in zip (X,Y):
f.write("%e %e \n"%(x,y))
f=None
def main():
derived_pars=Parameters.get_derived_pars()
print " RESOLUTION : " , derived_pars.step
xs=numpy.linspace(-Parameters.diam/2.0,Parameters.diam/2.0,Parameters.N, endpoint=False).astype(numpy.float32)
tmp = build_FZ(Parameters, derived_pars , xs )
FZP_2D = tmp["FZP_2D"]
CS_2D = tmp["CS_2D" ]
OSA_2D = tmp["OSA_2D"]
tmp=None
if(0) :
scrivi_linea(xs, FZP_2D, "linea_fz", Parameters.N//2,1)
trasmittance_FZP_CS = numpy.exp(-1j*FZP_2D*Parameters.phase) * CS_2D #
N_V_1 = int(Parameters.N/2-1-Parameters.slit_v/2.0/derived_pars.step)
N_V_2 = int(Parameters.N/2-1+Parameters.slit_v/2.0/derived_pars.step)
N_H_1 = Parameters.N/2-1+(Parameters.off_set-Parameters.slit_h/2.0)/derived_pars.step
N_H_2 = Parameters.N/2-1+(Parameters.off_set+Parameters.slit_h/2.0)/derived_pars.step
FIELD =numpy.zeros((Parameters.N,Parameters.N),dtype=numpy.complex64) # field_in_slit
FIELD[ N_V_1:N_V_2 , N_H_1:N_H_2 ] = 1.0
print " #### Propagation from slit position to FZP ####"
qx=( fft.fftfreq( Parameters.N ) * (2*math.pi*Parameters.N/Parameters.diam) ) .astype(numpy.float32)
q_2D = (qx*qx)[:,None]+(qx*qx)[None,:]
FIELD=fft.fftn(FIELD)
proj=(numpy.exp((-1.0j)*Parameters.slit_far*q_2D/(2*derived_pars.k))) .astype(numpy.complex64) #Fresnel propagator (Bunk&Diaz)
FIELD=fft.ifftn(proj*FIELD)
print "Propagated Complex Field from Slit : DONE !!!"
FIELD=FIELD*trasmittance_FZP_CS
print "Initial Complex Field for Partial illumination : DONE !!!"
print " ### Propagation of complex field (+CS+SLIT) to OSA ###"
FIELD=fft.fftn(FIELD)
proj=(numpy.exp((-1j)*derived_pars.OSA_position*q_2D/(2*derived_pars.k))) .astype(numpy.complex64) #Fresnel propagator (Bunk&Diaz)
FIELD=fft.ifftn(proj*FIELD)*OSA_2D
#numpy.save("/data/id01/inhouse/fra_quiney/FZP_70nmOZW/field_osa_in_2D_res20nm.npy",field_OSA_in)
print "Complex Field @ OSA : DONE !!!"
FIELD=fft.fftn(FIELD)
proj=(numpy.exp((-1j)* (derived_pars.focus-derived_pars.OSA_position) *q_2D/(2*derived_pars.k))) .astype(numpy.complex64) #Fresnel propagator (Bunk&Diaz)
FIELD=fft.ifftn(proj*FIELD)
if Parameters.FILE_FOR_FOCAL_DEPTH_PROFILING is not None:
print "#### looking at the longitudinal planes xz; yz ####"
x1=numpy.argmin(abs(xs+Parameters.range_x_focal_dep/2))
x2=numpy.argmin(abs(xs-Parameters.range_x_focal_dep/2))
Nsel= x2-x1
Ncorrection = roundfft(Nsel)-Nsel
x2=x2+(Ncorrection-Ncorrection/2)
x1=x1-(Ncorrection/2)
Nsel=Nsel+Ncorrection
xsel=xs[x1:x2]
FIELD_focus=FIELD[x1:x2,x1:x2]
qx=( fft.fftfreq( Nsel ) * (2*math.pi/((xsel[1]-xsel[0]))) ).astype(numpy.float32)
q_2D = (qx*qx)[:,None]+(qx*qx)[None,:]
z=numpy.linspace(-Parameters.range_z_focal_dep/2,Parameters.range_z_focal_dep/2, Parameters.Nz_focal_dep)
FFT_FIELD_focus= fft.fftn(FIELD_focus)
map_complex=numpy.zeros((Parameters.Nz_focal_dep,Nsel,Nsel),dtype=numpy.complex64)
for i in xrange(0,len(z)-1):
map_complex[i,:,:]= fft.ifftn(numpy.exp( -1j*z[i]*q_2D/(2*derived_pars.k)) * FFT_FIELD_focus)
print " writing results for focal depth profiling to " , Parameters.FILE_FOR_FOCAL_DEPTH_PROFILING+".h5"
h5=h5py.File( Parameters.FILE_FOR_FOCAL_DEPTH_PROFILING+".h5",'w')
h5["map_squareamplitude"]= map_complex.real*map_complex.real + map_complex.imag*map_complex.imag
h5["map_phase"] = numpy.arctan2(map_complex.imag,map_complex.real)
h5["map_real"] =map_complex.real
h5["map_imag"] =map_complex.imag
h5["z"]= z
h5["x"]= xsel
h5["y"]= xsel
h5.flush()
h5.close()
h5=None
####################################################################################
x1=numpy.argmin(numpy.abs(xs+derived_pars.range_x_focus2dect/2 ))
x2=numpy.argmin(numpy.abs(xs-derived_pars.range_x_focus2dect/2 ))+1
print "-------------------------- FOCUS 2 DECT --------------------- "
print " range_x_focus2dect = " , derived_pars.range_x_focus2dect
print " x1,x2 ", x1,x2
xsel=xs[x1:x2]
print " effective range_x_focus2dect (xsel[x2-1]-xsel[x1]) " , xs[x2-1]-xs[x1]
Nsel = len(xsel)
print " ##### calculating pixel size at the detector #####"
q_det=( fft.fftfreq( Nsel ) * (2*math.pi/((xs[1]-xs[0])))) .astype(numpy.float32)
angle_det=numpy.arcsin(derived_pars.lambd*q_det/2.0/math.pi).astype(numpy.float32)
x_det=Parameters.z_fd*numpy.tan(angle_det)
pxl=abs(x_det.max()-x_det.min())/(len(xsel)-1)
FIELD_select=FIELD[x1:x2,x1:x2]
print " #### propagator according to Fresnel free space approximation #### "
prop=(xsel*xsel)[:,None]+xsel*xsel
print " #### calculating compplex field at the detector #### "
B_fd=numpy.exp((+1j)*math.pi*prop/derived_pars.lambd/Parameters.z_fd).astype(numpy.complex64)
FIELD_select=fft.fftshift(fft.fftn( fft.ifftshift(B_fd*FIELD_select) ))
## this just changes an overall phase Q_d=(-1j)*FFT_fd*exp(1j*2*pi*z_fd/lambd)
#FZP_2D=numpy.load("/data/id01/inhouse/fra_quiney/FZP_70nmOZW/FZP_2D_res20nm.npy")
# phase0=0 ### Parameter phase slit ###
if Parameters.FILE_FOR_DETECTOR is not None :
print " writing results to " , Parameters.FILE_FOR_DETECTOR+".h5"
h5=h5py.File( Parameters.FILE_FOR_DETECTOR+".h5",'w')
h5["map_real"] =FIELD_select.real
h5["map_phase"] = numpy.arctan2(FIELD_select.imag,FIELD_select.real)
h5["map_imag"] =FIELD_select.imag
h5["map_squareamplitude"]= FIELD_select.real*FIELD_select.real + FIELD_select.imag*FIELD_select.imag
h5.flush()
h5.close()
h5=None
if(sys.argv[0][-12:]!="sphinx-build"):
main()