DUV Tutorial Series

# Install DUV simulation package
pip install duv-simulator

# Import required modules
import duv_simulator as duv
import numpy as np
import matplotlib.pyplot as plt

# Initialize simulator
simulator = duv.DUVSimulator()

# Basic Monte Carlo simulation
mc_results = simulator.run_monte_carlo(
    num_particles=100000,
    wavelength=193,  # nm
    na=0.85,
    feature_size=50.0,  # nm
    dose=30.0  # mJ/cm²
)

# Display results
print(f"Total energy deposited: {mc_results.statistics.total_energy:.2e} J")
print(f"Peak intensity: {mc_results.statistics.peak_intensity:.2e} W/cm²")
print(f"Simulation time: {mc_results.statistics.simulation_time:.2f} s")

# Create energy deposition heatmap
plt.figure(figsize=(10, 8))
plt.imshow(mc_results.energy_deposition, cmap='viridis')
plt.colorbar(label='Energy Density')
plt.title('Monte Carlo Energy Deposition Profile')
plt.xlabel('X Position (nm)')
plt.ylabel('Y Position (nm)')
plt.show()

# Plot particle distribution
plt.figure(figsize=(10, 6))
plt.scatter(mc_results.particle_positions[:, 0], 
           mc_results.particle_positions[:, 1],
           c=mc_results.particle_energies, 
           cmap='viridis', s=1)
plt.colorbar(label='Particle Energy')
plt.title('Particle Distribution')
plt.xlabel('X Position (nm)')
plt.ylabel('Y Position (nm)')
plt.show()

# Statistical analysis
statistics = mc_results.get_statistics()
print(f"Mean energy per particle: {statistics.mean_energy:.2e} J")
print(f"Energy standard deviation: {statistics.std_energy:.2e} J")
print(f"Energy range: {statistics.min_energy:.2e} - {statistics.max_energy:.2e} J")

# Parameter optimization
def optimize_simulation():
    # Test different particle counts
    particle_counts = [10000, 50000, 100000, 500000, 1000000]
    results = []
    
    for num_particles in particle_counts:
        result = simulator.run_monte_carlo(
            num_particles=num_particles,
            wavelength=193,
            na=0.85,
            feature_size=50.0,
            dose=30.0,
            simulation_speed="fast"
        )
        results.append({
            'particles': num_particles,
            'time': result.statistics.simulation_time,
            'energy': result.statistics.total_energy,
            'peak': result.statistics.peak_intensity
        })
    
    return results

# Run optimization
opt_results = optimize_simulation()

# Plot optimization results
plt.figure(figsize=(12, 4))

plt.subplot(1, 3, 1)
plt.plot([r['particles'] for r in opt_results], 
         [r['time'] for r in opt_results], 'o-')
plt.xlabel('Number of Particles')
plt.ylabel('Simulation Time (s)')
plt.title('Simulation Time vs Particles')

plt.subplot(1, 3, 2)
plt.plot([r['particles'] for r in opt_results], 
         [r['energy'] for r in opt_results], 'o-')
plt.xlabel('Number of Particles')
plt.ylabel('Total Energy (J)')
plt.title('Energy vs Particles')

plt.subplot(1, 3, 3)
plt.plot([r['particles'] for r in opt_results], 
         [r['peak'] for r in opt_results], 'o-')
plt.xlabel('Number of Particles')
plt.ylabel('Peak Intensity (W/cm²)')
plt.title('Peak Intensity vs Particles')

plt.tight_layout()
plt.show()

# Calculate double Gaussian PSF
psf_results = simulator.calculate_double_gaussian(
    sigma1=30.0,      # nm - narrow Gaussian
    sigma2=100.0,     # nm - wide Gaussian
    amplitude1=0.8,   # amplitude of narrow Gaussian
    amplitude2=0.2,   # amplitude of wide Gaussian
    grid_size=128,    # grid resolution
    feature_size=200.0  # nm - simulation area
)

# Display PSF metrics
print(f"FWHM1: {psf_results.metrics.fwhm1:.1f} nm")
print(f"FWHM2: {psf_results.metrics.fwhm2:.1f} nm")
print(f"Peak intensity: {psf_results.metrics.peak_intensity:.3f}")
print(f"Total energy: {psf_results.metrics.total_energy:.3f}")

# Visualize 2D PSF
plt.figure(figsize=(12, 5))

plt.subplot(1, 2, 1)
plt.imshow(psf_results.psf_2d, cmap='viridis')
plt.colorbar(label='PSF Value')
plt.title('Double Gaussian PSF (2D)')
plt.xlabel('X Position (nm)')
plt.ylabel('Y Position (nm)')

# Visualize 1D cross-section
plt.subplot(1, 2, 2)
x_positions = psf_results.psf_1d['x']
psf_values = psf_results.psf_1d['y']
plt.plot(x_positions, psf_values, 'b-', linewidth=2, label='Total PSF')

# Plot individual components
narrow_component = psf_results.narrow_component
wide_component = psf_results.wide_component
plt.plot(x_positions, narrow_component, 'r--', label='Narrow Gaussian')
plt.plot(x_positions, wide_component, 'g--', label='Wide Gaussian')

plt.xlabel('Position (nm)')
plt.ylabel('PSF Value')
plt.title('PSF Components (1D Cross-Section)')
plt.legend()
plt.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# Create test pattern
pattern = duv.create_line_pattern(
    width=50.0,  # nm
    length=200.0,  # nm
    grid_size=128
)

# Calculate aerial image
aerial_image = simulator.calculate_aerial_image(
    pattern=pattern,
    psf=psf_results.psf_2d
)

# Visualize results
plt.figure(figsize=(15, 5))

plt.subplot(1, 3, 1)
plt.imshow(pattern, cmap='gray')
plt.title('Test Pattern')
plt.xlabel('X Position (nm)')
plt.ylabel('Y Position (nm)')

plt.subplot(1, 3, 2)
plt.imshow(psf_results.psf_2d, cmap='viridis')
plt.title('Double Gaussian PSF')
plt.xlabel('X Position (nm)')
plt.ylabel('Y Position (nm)')

plt.subplot(1, 3, 3)
plt.imshow(aerial_image, cmap='viridis')
plt.title('Aerial Image')
plt.xlabel('X Position (nm)')
plt.ylabel('Y Position (nm)')

plt.tight_layout()
plt.show()