11"""A custom XY model for a supported DSPC bilayer."""
22
3- import math
3+ from math import sqrt
44
55import numpy as np
6+ from scipy .special import erf
67
78
89def custom_XY_DSPC (params , bulk_in , bulk_out , contrast ):
@@ -51,10 +52,10 @@ def custom_XY_DSPC(params, bulk_in, bulk_out, contrast):
5152 z = np .arange (0 , 141 )
5253
5354 # Make our Silicon substrate
54- vfSilicon , siSurf = layer (z , - 25 , 50 , 1 , subRough , subRough )
55+ vfSilicon , siSurf = make_layer (z , - 25 , 50 , 1 , subRough , subRough )
5556
5657 # Add the Oxide
57- vfOxide , oxSurface = layer (z , siSurf , oxideThick , 1 , subRough , subRough )
58+ vfOxide , oxSurface = make_layer (z , siSurf , oxideThick , 1 , subRough , subRough )
5859
5960 # We fill in the water at the end, but there may be a hydration layer between the bilayer and the oxide,
6061 # so we start the bilayer stack an appropriate distance away
@@ -65,15 +66,15 @@ def custom_XY_DSPC(params, bulk_in, bulk_out, contrast):
6566 headThick = vHead / lipidAPM
6667
6768 # ... and make a box for the volume fraction (1 for now, we correct for coverage later)
68- vfHeadL , headLSurface = layer (z , watSurface , headThick , 1 , bilayerRough , bilayerRough )
69+ vfHeadL , headLSurface = make_layer (z , watSurface , headThick , 1 , bilayerRough , bilayerRough )
6970
7071 # ... also do the same for the tails
7172 # We'll make both together, so the thickness will be twice the volume
7273 tailsThick = (2 * vTail ) / lipidAPM
73- vfTails , tailsSurf = layer (z , headLSurface , tailsThick , 1 , bilayerRough , bilayerRough )
74+ vfTails , tailsSurf = make_layer (z , headLSurface , tailsThick , 1 , bilayerRough , bilayerRough )
7475
7576 # Finally the upper head ...
76- vfHeadR , headSurface = layer (z , tailsSurf , headThick , 1 , bilayerRough , bilayerRough )
77+ vfHeadR , headSurface = make_layer (z , tailsSurf , headThick , 1 , bilayerRough , bilayerRough )
7778
7879 # Making the model
7980 # We've created the volume fraction profiles corresponding to each of the groups.
@@ -114,12 +115,12 @@ def custom_XY_DSPC(params, bulk_in, bulk_out, contrast):
114115 totSLD = sldSilicon + sldOxide + sldHeadL + sldTails + sldHeadR + sldWat
115116
116117 # Make the SLD array for output
117- SLD = [[ a , b ] for ( a , b ) in zip (z , totSLD )]
118+ SLD = np . column_stack ( (z , totSLD ))
118119
119120 return SLD , subRough
120121
121122
122- def layer (z , prevLaySurf , thickness , height , Sigma_L , Sigma_R ):
123+ def make_layer (z , prevLaySurf , thickness , height , Sigma_L , Sigma_R ):
123124 """Produce a layer, with a defined thickness, height and roughness.
124125
125126 Each side of the layer has its own roughness value.
@@ -129,12 +130,9 @@ def layer(z, prevLaySurf, thickness, height, Sigma_L, Sigma_R):
129130 right = prevLaySurf + thickness
130131
131132 # Make our heaviside
132- a = ( z - left ) / (( 2 ** 0.5 ) * Sigma_L )
133- b = ( z - right ) / (( 2 ** 0.5 ) * Sigma_R )
133+ erf_left = erf (( z - left ) / (sqrt ( 2 ) * Sigma_L ) )
134+ erf_right = erf (( z - right ) / (sqrt ( 2 ) * Sigma_R ) )
134135
135- erf_a = np .array ([math .erf (value ) for value in a ])
136- erf_b = np .array ([math .erf (value ) for value in b ])
137-
138- VF = np .array ((height / 2 ) * (erf_a - erf_b ))
136+ VF = np .array ((0.5 * height ) * (erf_left - erf_right ))
139137
140138 return VF , right
0 commit comments