# -*- coding: utf-8 -*-
"""
"""
import numpy
import math
import cPickle
from numpy import fft
import h5py
import sys
import h5py
def writeh5R(map_complex, nome):
h5=h5py.File( nome+".h5",'w')
h5["/map/map_real"] =map_complex
h5.flush()
h5.close()
h5=None
###### 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 :
""" TODO : document parameters """
E=8.0
diam=200*10**(-6)
delta_min=70*10**(-9)
N_pxl = 129
detector_pixel_size = 55.0e-6
pixels_per_degree = 300.0
supp_radius=302.0e-6
z_fd=None
###################################################################
datafile="data11Dec.npy"
previous_result_file = None
Beta=.85 ### HIO parameter always <1 ###
cycles_TOT=20
cycles_ER=500
cycles_HIO=1000
cycles_CF=500
thl_sw=8e-5 ### threshold for shrink-wrap
if(sys.argv[0][-12:]!="sphinx-build"):
exec(open(sys.argv[1],"r").read())
print "----- READAPTING N_pxl for fft "
print "..... OLD:",N_pxl,
# N_pxl=roundfft(N_pxl)
print "...... NEW:",N_pxl
allowed_pars = [ 'E', 'N_pxl', 'delta_min', 'diam',
"detector_pixel_size", "pixels_per_degree",
'datafile','previous_result_file', 'supp_radius',
'Beta','cycles_TOT','cycles_ER','cycles_HIO','cycles_CF','thl_sw',
'OSA_prox_factor', 'z_fd']
dic={}
exec(open(sys.argv[1],"r").read() ,globals(),dic)
extra_pars = set(dic.keys()) - set(allowed_pars+["allowed_pars"])
if extra_pars != set([]) :
print " PROBLEM: extra parameters have been set by input :", extra_pars
print " ------------------------------------------------------------- "
print " THE PARAMETERS I NEED ARE "
print allowed_pars
raise Exception,""
print " ########################################################################################## "
print " ### Now printing parameters. Stars means not set by input file (ERROR) ################# "
print " "
warn=0
for key in allowed_pars:
if key not in locals():
print "ERROR key " , key , " is missing "
raise Exception , " missing key "
if key in dic.keys():
print key," = " , locals()[key]
else:
print key," = " , locals()[key]," ******ERROR!!! *****"
warn=1
if warn:
print " ------------------------------------------------------------- "
print " THE PARAMETERS I NEED ARE "
print allowed_pars
print ""
print " use None for parameters like previous_result_file if they are not ready"
print ""
raise Exception,""
print " \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\n"
@staticmethod
[docs] def get_derived_pars():
""" TODO : document derived parameters """
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.r_max=self.diam/2
if self.z_fd is None:
del_theta_pxl=(1.0/self.pixels_per_degree ) * numpy.pi/180
res.z_fd =(self.pixels_per_degree)*self.detector_pixel_size*180.0/numpy.pi
res.size_focal_pixel=res.lambd/self.N_pxl/(del_theta_pxl) # size of a pixel at the focal plane
res.range_focus=self.N_pxl*res.size_focal_pixel
else :
res.z_fd = self.z_fd
self.pixels_per_degree = res.z_fd/( self.detector_pixel_size*180.0/numpy.pi )
del_theta_pxl=(1.0/self.pixels_per_degree ) * numpy.pi/180
res.size_focal_pixel=res.lambd/self.N_pxl/(del_theta_pxl)
res.range_focus=self.N_pxl*res.size_focal_pixel
span_q_det = 2*math.pi/res.size_focal_pixel
q_det = span_q_det * fft.fftfreq(self.N_pxl )
angle_det = numpy.arcsin(res.lambd*q_det/2/math.pi)
x_det = res.z_fd*numpy.tan(angle_det)
res.pxl_det= abs(x_det.max()-x_det.min())/(len(x_det)-1)
x_det = res.focus*numpy.tan(angle_det)
res.pxl_FZP= abs(x_det.max()-x_det.min())/(len(x_det)-1)
return res
def create_arrays_space_as_dictionary():
self=Parameters()
# Y0,X0=self.data_corner_YX
# datas=numpy.array(numpy.load(self.datafile) [Y0:Y0+self.N_pxl,X0:X0+self.N_pxl]).astype(numpy.float32)
datas=numpy.array(numpy.load(self.datafile)) .astype(numpy.float32)
if self.previous_result_file is not None:
solution=(numpy.load(self.previous_result_file)).astype(numpy.complex64)
else:
solution=numpy.zeros([self.N_pxl,self.N_pxl],numpy.complex64)
if not (solution.shape==datas.shape):
print " solution.shape==datas.shape, now padding to optimal N_pxl "
print solution.shape, datas.shape
datas_orig=datas
d0,d1 = datas_orig.shape
datas = numpy.zeros([self.N_pxl,self.N_pxl],numpy.float32)
datas[:d0,:d1] = datas_orig
datas[:d0,d1:].T[:,:]= datas_orig[:d0,d1-1]
datas[d0:,:d1]= datas_orig[d0-1,:d1]
datas[d0:,d1:]= datas_orig[d0-1,d1-1]
return {"datas":datas,"solution":solution}