refnxtoolbox/MetaStructure.py at master · igresh/refnxtoolbox · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157

import numpy as np

from scipy.interpolate import PchipInterpolator as Pchip
from scipy.integrate import simps

from refnx.reflect import Structure, Component, SLD, Slab
from refnx.analysis import Parameters, Parameter, possibly_create_parameter

import warnings

class metaComponent (Component):
    def __init__ (self, ComponentOne, ComponentTwo, microslab_max_thickness=1,
                  left_roughness=0, name='MixedLayer', ComponentOneSLD=None, ComponentTwoSLD=None):
        self.ComponentOne = ComponentOne
        self.ComponentTwo = ComponentTwo
        self.ComponentOneSLD = ComponentOneSLD
        self.ComponentTwoSLD = ComponentTwoSLD

        self._interfaces = None

        self.microslab_max_thickness = microslab_max_thickness

        self.left_roughness = left_roughness

        self.name = name


    def __call__ (self, z=None):
        zed1, vfp1 = self.get_vfp(self.ComponentOne.slabs())
        zed2, vfp2 = self.get_vfp(self.ComponentTwo.slabs())

        zed, vfp1, vfp2 = self.establish_basis (zed1, vfp1, zed2, vfp2)

        if np.any(z):
            vfp1 = np.interp(z, zed, vfp1)
            vfp2 = np.interp(z, zed, vfp2)
        else:
            z = zed

        return z, vfp1, vfp2


    def combined_sld(self, z=None):
        """
        Really just for testing


        """
        z, vfp1, vfp2 = self(z)
        total_vfp = vfp1 + vfp2

        if np.any(total_vfp>1.001):
            print (f'your total vf is greater than one: {np.max(total_vfp)}')


        if self.ComponentOneSLD is None:
            sld1 = self.ComponentOne.sld.real.value
        else:
            sld1 = self.ComponentOneSLD.real.value

        if self.ComponentTwoSLD is None:
            sld2 = self.ComponentTwo.sld.real.value
        else:
            sld2 = self.ComponentTwoSLD.real.value

        norm_vfp1 = vfp1/total_vfp
        norm_vfp2 = vfp2/total_vfp

        sld12 = sld1*norm_vfp1 + sld2*norm_vfp2


        return z, sld12, total_vfp


    def slabs(self, structure=None):

        z_init, vfp1_init, vfp2_init = self()
        slab_extent = z_init[-1]

        num_slabs = np.ceil(float(slab_extent) / self.microslab_max_thickness)
        slab_thick = float(slab_extent / num_slabs)
        slabs = np.zeros((int(num_slabs), 5))
        slabs[:, 0] = slab_thick

        # give last slab a miniscule roughness so it doesn't get contracted
        slabs[-1:, 3] = 0.5


        dist = np.cumsum(slabs[..., 0]) - 0.5 * slab_thick

        z, comb_sld, total_vfp = self.combined_sld(z=dist)

        slabs[:, 1] = comb_sld
        slabs[:, 2] = 0
        slabs[:, 4] = 1-total_vfp

        slabs[0,3] =  self.left_roughness

        return slabs


    @property
    def parameters(self):
        p = Parameters(name=self.name)
        p.extend(self.ComponentOne.parameters)
        p.extend(self.ComponentTwo.parameters)
        return p


    def logp(self):
        return 0


    def get_vfp(self, slabs):
        slabs = np.atleast_2d(slabs)
        s = Structure()
        s |= SLD(0)

        m = SLD(1.)

        for i, slab in enumerate(slabs):
            layer = m(slab[0], slab[3])
            layer.vfsolv.value = slab[4]
            s |= layer

        final_layer=SLD(0)(0)
        final_layer.vfsolv.setp(1)
        s |= final_layer

        s.solvent = SLD(0)


        total_thickness = np.sum(s.slabs()[:, 0])

        # Create a z axis for the volume fraction profile, with a
        # value every angstrom
        zed = np.linspace(0, total_thickness, int(total_thickness+1))

        z, s = s.sld_profile(z=zed)

        return z, s

    def establish_basis(self, zed1, vfp1, zed2, vfp2):

        high_zed = np.max([zed1[-1], zed2[-1]])

        new_zed = np.linspace(0, int(high_zed), int(high_zed+1))

        new_vfp1 = np.interp(new_zed, zed1, vfp1, right=0)
        new_vfp2 = np.interp(new_zed, zed2, vfp2, right=0)

        return new_zed, new_vfp1, new_vfp2