<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
Step1:  Residual Stress Intensity Factor - $K_{res}$
Background
The stress intensity factor for symettrical loading at the crack tip (x=a) is given by
$$K_I = \int^a_0 ! \sigma (x) \, h(a, x) \, \mathrm{d}x $$
where $\sigma(x)$ is the stress in an uncracked body and h(a,x) the weight function that is independent of $\sigma(x)$. To calculate $K_{res}$ we therefore need accurate expressions for $\sigma_{res}(x)$ and $h_{CT}(a, x)$.
Weight Function - h(a, x)
A number of researchers have come up with weight function approximations for different geometries (see weight function notebook). In this instace we are primarily interested in the compact-tension (CT) geometry. One such approximation was provided by Fett and Munz (1997) who showed that
$$ h(a, x) = \frac{1}{\sqrt{2 \pi a}} 
\frac{1}{\sqrt{1 - \frac{x}{a}}}
\left[1 + \sum^\mu_v
\frac{A_{v, \mu}\left( \frac{a}{W}\right)^2}
{\left(1- \frac{a}{W}\right)^{\frac{3}{2}}}
\left(1 - \frac{x}{a}
\right)^{v+1}
\right]$$
where
Step2:  We can evaluate that weight function at different a/W combinations
Step3:  Residual Stress - $\sigma_{res}$
We need to extract data and calculate strains and stresses from an experiment carried out at ID11. In particular, we'll interrogate data taken prior to pre-cracking a bainitic compact-tension sample subject to a compressive pre-strain. 
Reload Diffraction Data
This data initially came as a folder containing 700+ images (Debye-Scherrer rings) and a SPEC metadata file. Positional information was scraped from the SPEC file and azimuthal integration was completed on the associated image files using pyXe and the underlying pyFAI functions. The 1D line profiles (as a function of phi) were saved to a HDF5 file and will be reloaded here.
Prior to testing a far field scan was completed with the aim of extracting valid d0/q0 values.
Step4:  Calculate Strain (apply q0/d0)
Static q0/d0
Step5:  Apply q0 scan
Should be able to simply pass the q0 points to the test data and the phi normalised strain should be accurately calculated. The best way to test the resultant values is to plot the final strain fit, which should vary sinosoidal according to the strain transformation equation
Step6:  This is better but surprisingly still not fantastically good - it raises questions over whether the d0/q0 strain free acquisition was as good as it needed to be. We should try to do better
Step7:  Both are much better - at this stage is seems more reasonable to use the far field values (although the fit is slightly worse).
Sample Geometry
Before extracting the 1D line profile we need to align the 2D strain map, such that the x=0 reflects the same position as in the weight function calculation. To do this we first need to know the geometry of the sample
Step8:  Stress Profile Extraction
We need to extract the stress profile, rather than the strain profile. To do this we need to apply the materials constants - i.e. poisson's ratio and Young's modulus.
Step9:  Clean Data
The data at approximately x = 10mm looks suspect. From the 2D map we can also see that this region also sees a peak compressive strain. This looks to be some kind of partial volume/edge effect, with the data being nonsensical (and effectively equal to 0). We can attempt to test this hypothesis by integrating the line profile. Given that the stress returns to zero in the far field, the integrated stress should also be zero. If the peak at 10mm is incorrect then this will raise the integrated stress above 0.
Step10:  As expected, the integrated stress is closer to zero when the start point is taken to be the second peak. The integrated stress is 0 when x = 11.4mm, although for now we will leave x0 = 11mm. If we consider the data point to be invalid (in air) then lets make z = 0MPa when x < 11mm.
Step11:  Calculate $K_{res}$
We're now finally in a position where we can calculate $K_{res}$ as a function of position along the residual stress profile.
Step12:  Save Data
<ASSISTANT_TASK:>
Python Code:
def weight_CT_FettMunz(a, W, dx=0.005):
    Weight function estimation as per Fett and Munz (1997).
    
    Args:
        a (float): Crack length
        W (float): Width of CT sample
        dx (float): Point spacing for calculations
        
    Returns:
        tuple: x, weight - position and weighting of point
    
    
    A_vmu = np.array([[2.673, -8.604, 20.621, -14.635, 0.477],
                      [-3.557, 24.973, -53.398, 50.707, -11.837],
                      [1.230, -8.411, 16.957, -12.157, -0.940],
                      [-0.157, 0.954, -1.284, -0.393, 1.655]])
    
    # Remove if x[-1] if it equals a (due to zero div error)
    x = np.arange(0, a, dx)
    x = x if x[-1] != a else x[:-1]
    rho = x / a
    aW = a/W
    h1 = (2 / (np.pi * a))**0.5 * (1 / (1 - rho))**0.5
    h2 = 1 + x * 0
    for v, mu in np.ndindex(4,5):
        A = A_vmu[v, mu]
        h2 += A * aW**mu * (1 - rho)**(v + 1) / (1 - aW)**1.5
    return x, h1 * h2
Explanation: Residual Stress Intensity Factor - $K_{res}$
Background
The stress intensity factor for symettrical loading at the crack tip (x=a) is given by
$$K_I = \int^a_0 ! \sigma (x) \, h(a, x) \, \mathrm{d}x $$
where $\sigma(x)$ is the stress in an uncracked body and h(a,x) the weight function that is independent of $\sigma(x)$. To calculate $K_{res}$ we therefore need accurate expressions for $\sigma_{res}(x)$ and $h_{CT}(a, x)$.
Weight Function - h(a, x)
A number of researchers have come up with weight function approximations for different geometries (see weight function notebook). In this instace we are primarily interested in the compact-tension (CT) geometry. One such approximation was provided by Fett and Munz (1997) who showed that
$$ h(a, x) = \frac{1}{\sqrt{2 \pi a}} 
\frac{1}{\sqrt{1 - \frac{x}{a}}}
\left[1 + \sum^\mu_v
\frac{A_{v, \mu}\left( \frac{a}{W}\right)^2}
{\left(1- \frac{a}{W}\right)^{\frac{3}{2}}}
\left(1 - \frac{x}{a}
\right)^{v+1}
\right]$$
where:
<img src="Images/CT_FettMunz.png" width="700"/>
End of explanation
# Calculate the weight function at a/W = 0.3, 0.5, 0.7
plt.figure(figsize=(7, 5))
for linspec, a in zip(['k:', 'k--', 'k'], 
                      [0.3, 0.5, 0.7]): 
    x, h = weight_CT_FettMunz(a, 1, dx=0.0025)
    plt.plot(x/a, h, linspec, label='a/W = {}'.format(a))
plt.legend()
plt.xlabel('Relative Position, x/a')
plt.ylabel('Weight function, h(a, x)')
plt.show()
Explanation: We can evaluate that weight function at different a/W combinations:
End of explanation
base_folder = os.path.expanduser('~/Dropbox/Synchrotron Data/id11')
# First extract the q0 values
q0_loc = r'd0_2_/d0_2__pyxe.h5'
file = os.path.join(base_folder, q0_loc)
q0 = PeakAnalysis(file)
q0.peak_fit(3.1, 0.5)
# Then open the real data
data_loc = r'thick4_coarse_0b/thick4_coarse_0b_pyxe.h5'
file = os.path.join(base_folder, data_loc)
data = PeakAnalysis(file)
# data.peak_fit(3.1, 0.5)
Explanation: Residual Stress - $\sigma_{res}$
We need to extract data and calculate strains and stresses from an experiment carried out at ID11. In particular, we'll interrogate data taken prior to pre-cracking a bainitic compact-tension sample subject to a compressive pre-strain. 
Reload Diffraction Data
This data initially came as a folder containing 700+ images (Debye-Scherrer rings) and a SPEC metadata file. Positional information was scraped from the SPEC file and azimuthal integration was completed on the associated image files using pyXe and the underlying pyFAI functions. The 1D line profiles (as a function of phi) were saved to a HDF5 file and will be reloaded here.
Prior to testing a far field scan was completed with the aim of extracting valid d0/q0 values.
End of explanation
q0_mean = np.nanmean(q0.peaks)
print('q0 = {:.4f}'.format(q0_mean))
data.calculate_strain(q0_mean)
data.plot_strain_fit(pnt=100)
Explanation: Calculate Strain (apply q0/d0)
Static q0/d0
End of explanation
data.calculate_strain(q0)
data.plot_strain_fit(pnt=100)
Explanation: Apply q0 scan
Should be able to simply pass the q0 points to the test data and the phi normalised strain should be accurately calculated. The best way to test the resultant values is to plot the final strain fit, which should vary sinosoidal according to the strain transformation equation:
End of explanation
data.calculate_strain(np.nanmean(data.peaks, 0))
data.plot_strain_fit(pnt=200)
# Mean from final 3 rows
data.calculate_strain(np.nanmean(data.peaks[-63:], 0))
data.plot_strain_fit(pnt=200)
Explanation: This is better but surprisingly still not fantastically good - it raises questions over whether the d0/q0 strain free acquisition was as good as it needed to be. We should try to do better:
q0 vector
We'll take two approaches to this:
Mean q profile from all data points. If we have sampled enough material (infront and behind crack) then the total strain will sum to zero and the average profile should be equivalent to the unstrained profile.
Mean q profile from further field values. Use points at the furthest extent of the profile. If they are sufficiently away from the crack tip/strain field then these can be considered to be far-field strain free measurements.
End of explanation
# Rotate and flip
data.swapaxes(0, 1)
data.flipaxis(1)
data.centre((-3.575-11, -10.7))
data.plot_slice('strain', phi=np.pi/2, figsize=(7,7))
Explanation: Both are much better - at this stage is seems more reasonable to use the far field values (although the fit is slightly worse).
Sample Geometry
Before extracting the 1D line profile we need to align the 2D strain map, such that the x=0 reflects the same position as in the weight function calculation. To do this we first need to know the geometry of the sample:
<img src="Images/sample.png" alt="Drawing" width="400"/>
In addition to these dimensions, it is important to note that that a small notch was then added to the pre-strained sample to aid crack initiation. This was approximately 1mm in length. Taking all of this into account, the notch tip is loacted at a length, x of 11mm (main notch 10mm in length). 
Data Transformation
With the geometry/notch in mind, lets plot the transformed (rotated/flipped/recentred) 2D strain map.
End of explanation
data.material_parameters(E=200e9, v=0.3)
plt.figure(figsize=(6, 5))
x, y, d, sigma = data.extract_line(data='stress', phi=np.pi/2, pnt=(11, 0))
sigma /= 1e6
x /= 1e3
plt.plot(x * 1e3, sigma)
x0 = [i/1e3 for i in [10.05, 10.5, 11.0]]
for i in x0:
    plt.axvline(i * 1e3, ls='--', color='k')
Explanation: Stress Profile Extraction
We need to extract the stress profile, rather than the strain profile. To do this we need to apply the materials constants - i.e. poisson's ratio and Young's modulus.
End of explanation
# Integrate from first peak, mid-point trough, and second peak
int_stress = []
for i in x0:
    valid = np.logical_and(np.isfinite(sigma), x > i) 
    int_stress.append(np.trapz(sigma[valid], x[valid]))
plt.figure(figsize=(7, 5))
plt.plot([i * 1e3 for i in x0], int_stress, 'k.-')
plt.xlabel('Integration start point, x_0 (mm)')
plt.ylabel('Stress sum (MPa)')
Explanation: Clean Data
The data at approximately x = 10mm looks suspect. From the 2D map we can also see that this region also sees a peak compressive strain. This looks to be some kind of partial volume/edge effect, with the data being nonsensical (and effectively equal to 0). We can attempt to test this hypothesis by integrating the line profile. Given that the stress returns to zero in the far field, the integrated stress should also be zero. If the peak at 10mm is incorrect then this will raise the integrated stress above 0.
End of explanation
dx = 0.00051 / 1e3
x_new = np.concatenate((np.array([0]), x[np.isfinite(sigma)]))
sigma_new = np.concatenate((np.array([0]), sigma[np.isfinite(sigma)]))
sigma_new[x_new < 11/1e3] = 0
f = interp1d(x_new, sigma_new, 1)
# Lets check that interpolation with spacings defined by weight function calcs
x_h = np.arange(0, x_new.max(), dx)
plt.figure(figsize=(7, 5))
plt.plot(x_new * 1e3, sigma_new, 'r')
sig_b = f(x_h)
plt.plot(x_h * 1e3, sig_b, 'k--')
plt.ylabel(r'Stress (MPa)')
plt.xlabel(r'Position (mm)')
Explanation: As expected, the integrated stress is closer to zero when the start point is taken to be the second peak. The integrated stress is 0 when x = 11.4mm, although for now we will leave x0 = 11mm. If we consider the data point to be invalid (in air) then lets make z = 0MPa when x < 11mm.
End of explanation
a = np.linspace(0, x_new.max(), 201)
K = np.nan * np.zeros_like(a)
for i, ai in enumerate(a):
    try:
        x_h, h = weight_CT_FettMunz(ai, 50/1e3, dx=dx)
        sigma_x = f(x_h)
        K[i] = np.trapz(h * sigma_x, x_h)
    except IndexError:
        pass
fig = plt.figure(figsize=(7, 5))
plt.ylabel(r'$\mathrm{K_{res}}$ $\mathrm{(MPam^{1/2})}$', size=14)
plt.xlabel(r'Position (mm)', size=14)
plt.plot(a * 1e3, K, 'k')
plt.xlim([10, 27])
Explanation: Calculate $K_{res}$
We're now finally in a position where we can calculate $K_{res}$ as a function of position along the residual stress profile.
End of explanation
f_end = r'~/Dropbox/Papers, Proposals and Presentations/Papers/ECF21/Figures'
path = os.path.join(os.path.expanduser(f_end), r'k_res_tens.csv')
print(path)
sig_a = f(a)
np.savetxt(path, (a, sig_a, K), delimiter=',')
len(K)
Explanation: Save Data
End of explanation
<END_TASK>