0

I am trying to write a function to perform a Simpson’s rule integration of the 1-dimensional Fresnel integral. And am as new as can be to coding. I am struggling to write a code the solves the integral and plots it for all possible values using arrays. More precisely, I am trying to figure out how to describe y where x is the argument and x' is in a loop of x' limits, which are the limits of the integral, as shown in the picture. Any further recommended changes to my code would be very helpful. The format of the question that was given

This is my current effort. Includes and indented block I cannot get rid of. 

import math 
import cmath
import numpy as np
import matplotlib.pyplot as plt

NumPoints = 100 #Number of areas to be calculated for integration           approximation
λ = 1*10**-6 #Wavelength 
z = 0.02 #Screen distance
k = (2*math.pi)/λ #wave number
j = j=cmath.sqrt(-1)

xmin = -0.005 #x = screen coordinate
xmax = +0.005 
xvals = np.linspace(xmin,xmax,NumPoints) # x intervals between limits

xdmin = 0.0 #xd = aperture coordinate 
xdmax = 2*10**-5
dxd = xdmax/NumPoints
xdvals = np.linspace(xdmin,xdmax,Numpoints) # xd intervals between limits
#All figures in SI

def simpson (x, xd ):

    for i in range(Numpoints):
    while (xdmin<xd<xdmax):

    yvals[i] = (k/2(math.pi)z)math.exp((1jk/2z)(xvals[i] - xd)**2)

    integral = np.sum(yvals[i])*dxd 

    return integral

print("The integral =", simpson, "for the given constants." )

plt.plot(xvals,yvals)

Thanks!

Mr Singh
  • 3,936
  • 5
  • 41
  • 60
user471264
  • 3
  • 4

1 Answers1

0

If it's not important to write your own Simpon's Rule, then you can use the import from scipy.integrate import simps. The following code is an attempt at your problem: 

import math
import cmath
from scipy.integrate import simps
import numpy as np

# Variables
wavelength = 1*10**(-6) # Wavelength
z = 0.02 # Screen distance
k = (2*math.pi) / wavelength # Wave number
n = 100
aperture_width = 2*10**(-5)

# Define the integrand
def f(a, b):
    integrand = cmath.exp(1j*k/(2*z)*(a - b)**2)
    return integrand

# Paramters defining the limits of integration and x'
x = np.linspace(-0.005, 0.005, n)
xd = np.linspace(0.0, aperture_width, n)

# Create a list of y values from the integrand function above
y = [0.0]*n
for i in range(1, n):
    y[i] = f(x[i], xd[i])

# Using Simpon's Rule, integrate y dx
print(k/(2*math.pi*z)*simps(y, x)) # Output is complex

Running this code produces the output (-26942.956248350412-72020.42454065284j). If this answer is too ugly, you can simplify the answer using scientific notation as follows: 

# Using Simpon's Rule, integrate y dx
answer = k/(2*math.pi*z)*simps(y, x)
simplified_answer = '{:.2e}'.format(answer)  # Round and put in scientific notation

print(simplified_answer)  # Output is complex

This results in a much prettier -2.69e+04-7.20e+04j.

Andy
  • 789
  • 8
  • 19