Source code for ParametersR_module

# -*- 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}