pyapi07_hexspheres_chi2.py

pospelov, 19 Dec 2018 08:23

Download (3.48 KB)

 
1
import numpy as np
2
import bornagain as ba
3
from bornagain import deg, angstrom, nm, millimeter
4
from matplotlib import pyplot as plt
5
from mpl_toolkits.mplot3d import Axes3D
6

    
7

    
8
def get_sample(params):
9
    radius = params["radius"]
10
    length = params["length"]
11

    
12
    # Defining Materials
13
    material_1 = ba.HomogeneousMaterial("Air", 0.0, 0.0)
14
    material_2 = ba.HomogeneousMaterial("CoFe2O4", 2.0433e-05, 1.5253e-06)
15
    material_3 = ba.HomogeneousMaterial("SiO2", 5.43852457e-06, 5.43741763e-08)
16
    material_4 = ba.HomogeneousMaterial("Si", 5.78164999998e-06, 1.02295e-07)
17

    
18
    # Defining Layers
19
    layer_1 = ba.Layer(material_1)
20
    layer_2 = ba.Layer(material_3, 60)
21
    layer_3 = ba.Layer(material_4)
22

    
23
    # Defining Form Factors
24
    formFactor_1 = ba.FormFactorTruncatedSphere(radius * nm, 8.5 * nm, 0.0 * nm)
25
    # formFactor_1 = ba.FormFactorFullSphere(radius)
26

    
27
    # Defining Particles
28
    particle_1 = ba.Particle(material_2, formFactor_1)
29

    
30
    # Defining Interference Functions
31
    interference_1 = ba.InterferenceFunction2DLattice(length, length, 120.0 * deg, 0.0 * deg)
32
    interference_1_pdf = ba.FTDecayFunction2DCauchy(300.0 * nm, 100.0 * nm, 0.0 * deg)
33
    interference_1.setDecayFunction(interference_1_pdf)
34

    
35
    # Defining Particle Layouts and adding Particles
36
    layout_1 = ba.ParticleLayout()
37
    layout_1.addParticle(particle_1, 1.0)
38
    layout_1.setInterferenceFunction(interference_1)
39
    layout_1.setTotalParticleSurfaceDensity(0.000998875898252)
40

    
41
    # Adding layouts to layers
42
    layer_1.addLayout(layout_1)
43

    
44
    # Defining Multilayers
45
    multiLayer_1 = ba.MultiLayer()
46
    multiLayer_1.addLayer(layer_1)
47
    multiLayer_1.addLayer(layer_2)
48
    multiLayer_1.addLayer(layer_3)
49
    return multiLayer_1
50

    
51

    
52
def get_simulation(params):
53
    simulation = ba.GISASSimulation()
54

    
55
    detector = ba.RectangularDetector(981, 168.732, 1043, 179.396)
56
    detector.setPerpendicularToDirectBeam(3532.0, 103.234, 60.062)
57
    simulation.setDetector(detector)
58

    
59
    simulation.setRegionOfInterest(87, 70, 119, 130)
60
    simulation.addMask(ba.Rectangle(101.7, 72, 104.7, 107), True)
61

    
62
    simulation.setDetectorResolutionFunction(ba.ResolutionFunction2DGaussian(0.1, 0.1))
63
    simulation.setBeamParameters(0.1 * nm, 0.2 * deg, 0.0 * deg)
64
    simulation.setBeamIntensity(1.0e+08)
65

    
66
    simulation.setSample(get_sample(params))
67
    return simulation
68

    
69

    
70
def plot_chi2():
71

    
72
    ref_data = np.loadtxt("experimental_data.txt.gz")
73
    a_radius = np.linspace(5.0, 15.0, num=200)
74
    a_length = np.linspace(20.0, 40.0, num=200)
75

    
76
    objective = ba.FitObjective()
77
    objective.addSimulationAndData(get_simulation, ref_data)
78

    
79
    chi2_data = np.zeros([len(a_radius), len(a_length)])
80
    for i_rad, radius in enumerate(a_radius):
81
        for i_len, length in enumerate(a_length):
82
            params = ba.Parameters()
83
            params.add("radius", value=radius)
84
            params.add("length", value=length)
85
            objective.evaluate(params)
86
            value = np.log10(objective.iterationInfo().chi2())
87
            chi2_data[i_rad][i_len] = value
88
            print("#{0} of {1} (radius={2:-6.2f}) (length={3:-6.2f}) (chi2={4:-6.2f})".format(i_rad, len(a_radius), radius, length, value))
89

    
90
    x, y = np.meshgrid(a_radius, a_length)
91

    
92
    fig = plt.figure(figsize=(12.80, 10.24))
93
    ax = fig.gca(projection='3d')
94
    ax.plot_surface(x, y, chi2_data, rstride=2, cstride=2, alpha=0.3)
95
    ax.set_xlabel('radius')
96
    ax.set_ylabel('length')
97
    ax.set_zlabel('log10(chi2)')
98
    ax.view_init(elev=20., azim=20)
99

    
100

    
101
if __name__ == '__main__':
102
    plot_chi2()
103
    plt.show()