import math\n",
"\n",
"import numpy as np\n",
"\n",
"\n",
"def domains_XY_model(params, bulk_in, bulk_out, contrast, domain):\n",
" # Split up the parameters for convenience\n",
" subRough = params[0]\n",
" oxideThick = params[1]\n",
" layerThick = params[2]\n",
" layerSLD = params[3]\n",
" layerRough = params[4]\n",
" domainSLD = params[5]\n",
"\n",
" # Make an array of z values for our model\n",
" z = np.arange(0, 141)\n",
"\n",
" # Make the volume fraction distribution for our Silicon substrate\n",
" [vfSilicon, siSurf] = makeLayer(z, -25, 50, 1, subRough, subRough)\n",
"\n",
" # ... and the Oxide ...\n",
" [vfOxide, oxSurface] = makeLayer(z, siSurf, oxideThick, 1, subRough, subRough)\n",
"\n",
" # ... and also our layer.\n",
" [vfLayer, laySurface] = makeLayer(z, oxSurface, layerThick, 1, subRough, layerRough)\n",
"\n",
" # Everything that is not already occupied will be filled will water\n",
" totalVF = vfSilicon + vfOxide + vfLayer\n",
" vfWater = 1 - totalVF\n",
"\n",
" # Now convert the Volume Fractions to SLDs\n",
" siSLD = vfSilicon * bulk_in\n",
" oxSLD = vfOxide * 3.41e-6\n",
"\n",
" # Layer SLD depends on whether we are calculating the domain or not\n",
" if domain == 0:\n",
" laySLD = vfLayer * layerSLD\n",
" else:\n",
" laySLD = vfLayer * domainSLD\n",
"\n",
" # ... and finally the water SLD.\n",
" waterSLD = vfWater * bulk_out[contrast]\n",
"\n",
" # Make the total SLD by just adding them all up\n",
" totalSLD = siSLD + oxSLD + laySLD + waterSLD\n",
"\n",
" # The output is just a [n x 2] array of z against SLD\n",
" SLD = np.column_stack([z, totalSLD])\n",
"\n",
" return SLD, subRough\n",
"\n",
"\n",
"def makeLayer(z, prevLaySurf, thickness, height, Sigma_L, Sigma_R):\n",
" """This produces a layer, with a defined thickness, height and roughness.\n",
" Each side of the layer has its own roughness value.\n",
" """\n",
" # Find the edges\n",
" left = prevLaySurf\n",
" right = prevLaySurf + thickness\n",
"\n",
" # Make our heaviside\n",
" a = (z - left) / ((2**0.5) * Sigma_L)\n",
" b = (z - right) / ((2**0.5) * Sigma_R)\n",
"\n",
" erf_a = np.array([math.erf(value) for value in a])\n",
" erf_b = np.array([math.erf(value) for value in b])\n",
"\n",
" VF = np.array((height / 2) * (erf_a - erf_b))\n",
"\n",
" return VF, right\n",
"