DUV Code Examples

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

class MonteCarloSimulator:
    def __init__(self, wavelength=193, na=0.85):
        self.wavelength = wavelength
        self.na = na
        self.particles = []
        
    def generate_particles(self, num_particles, feature_size, dose):
        """Generate particles with random positions and energies"""
        # Random positions within feature area
        x = np.random.uniform(-feature_size/2, feature_size/2, num_particles)
        y = np.random.uniform(-feature_size/2, feature_size/2, num_particles)
        
        # Random energies based on dose
        energies = np.random.exponential(dose * 1e-6, num_particles)
        
        # Random directions
        directions = np.random.uniform(0, 2*np.pi, num_particles)
        
        self.particles = np.column_stack([x, y, energies, directions])
        return self.particles
    
    def simulate_energy_deposition(self, grid_size=200):
        """Simulate energy deposition on grid"""
        energy_grid = np.zeros((grid_size, grid_size))
        pixel_size = 2.0 / grid_size  # 2 μm total area
        
        for particle in self.particles:
            x, y, energy, direction = particle
            
            # Convert to grid coordinates
            grid_x = int((x + 1.0) / pixel_size)
            grid_y = int((y + 1.0) / pixel_size)
            
            if 0 <= grid_x < grid_size and 0 <= grid_y < grid_size:
                energy_grid[grid_y, grid_x] += energy
        
        return energy_grid
    
    def calculate_statistics(self, energy_grid):
        """Calculate simulation statistics"""
        total_energy = np.sum(energy_grid)
        peak_intensity = np.max(energy_grid)
        mean_energy = np.mean(energy_grid[energy_grid > 0])
        std_energy = np.std(energy_grid[energy_grid > 0])
        
        return {
            'total_energy': total_energy,
            'peak_intensity': peak_intensity,
            'mean_energy': mean_energy,
            'std_energy': std_energy
        }
    
    def visualize_results(self, energy_grid):
        """Visualize simulation results"""
        fig, axes = plt.subplots(1, 2, figsize=(12, 5))
        
        # Energy deposition heatmap
        im1 = axes[0].imshow(energy_grid, cmap='viridis')
        axes[0].set_title('Energy Deposition Profile')
        axes[0].set_xlabel('X Position (μm)')
        axes[0].set_ylabel('Y Position (μm)')
        plt.colorbar(im1, ax=axes[0], label='Energy Density')
        
        # Particle distribution
        if len(self.particles) > 0:
            x_pos = self.particles[:, 0]
            y_pos = self.particles[:, 1]
            energies = self.particles[:, 2]
            
            scatter = axes[1].scatter(x_pos, y_pos, c=energies, 
                                    cmap='viridis', s=1, alpha=0.6)
            axes[1].set_title('Particle Distribution')
            axes[1].set_xlabel('X Position (μm)')
            axes[1].set_ylabel('Y Position (μm)')
            plt.colorbar(scatter, ax=axes[1], label='Particle Energy')
        
        plt.tight_layout()
        plt.show()

# Usage example
simulator = MonteCarloSimulator(wavelength=193, na=0.85)
particles = simulator.generate_particles(num_particles=100000, 
                                       feature_size=1.0, dose=30.0)
energy_grid = simulator.simulate_energy_deposition()
stats = simulator.calculate_statistics(energy_grid)
simulator.visualize_results(energy_grid)

print(f"Total energy: {stats['total_energy']:.2e} J")
print(f"Peak intensity: {stats['peak_intensity']:.2e} W/cm²")

import numpy as np
import matplotlib.pyplot as plt
from scipy.fft import fft2, ifft2, fftshift

class DoubleGaussianPSF:
    def __init__(self, grid_size=128, feature_size=2.0):
        self.grid_size = grid_size
        self.feature_size = feature_size
        self.pixel_size = feature_size / grid_size
        
    def calculate_psf(self, sigma1, sigma2, amplitude1, amplitude2):
        """Calculate double Gaussian PSF"""
        # Create coordinate grids
        x = np.linspace(-feature_size/2, feature_size/2, grid_size)
        y = np.linspace(-feature_size/2, feature_size/2, grid_size)
        X, Y = np.meshgrid(x, y)
        
        # Calculate distances from center
        R = np.sqrt(X**2 + Y**2)
        
        # Calculate double Gaussian
        gaussian1 = amplitude1 * np.exp(-R**2 / (2 * sigma1**2))
        gaussian2 = amplitude2 * np.exp(-R**2 / (2 * sigma2**2))
        psf = gaussian1 + gaussian2
        
        return psf, X, Y
    
    def calculate_aerial_image(self, pattern, psf):
        """Calculate aerial image using FFT convolution"""
        # FFT of pattern and PSF
        pattern_fft = fft2(pattern)
        psf_fft = fft2(psf)
        
        # Convolution in frequency domain
        aerial_fft = pattern_fft * psf_fft
        
        # Inverse FFT
        aerial_image = np.real(ifft2(aerial_fft))
        
        return aerial_image
    
    def calculate_metrics(self, psf):
        """Calculate PSF metrics"""
        # Find peak
        peak_intensity = np.max(psf)
        
        # Calculate FWHM
        center = psf.shape[0] // 2
        profile = psf[center, :]
        half_max = peak_intensity / 2
        
        # Find FWHM points
        above_half_max = profile > half_max
        if np.any(above_half_max):
            indices = np.where(above_half_max)[0]
            fwhm = (indices[-1] - indices[0]) * self.pixel_size
        else:
            fwhm = 0
        
        # Total energy
        total_energy = np.sum(psf) * self.pixel_size**2
        
        return {
            'peak_intensity': peak_intensity,
            'fwhm': fwhm,
            'total_energy': total_energy
        }
    
    def visualize_psf(self, psf, X, Y):
        """Visualize PSF results"""
        fig, axes = plt.subplots(1, 3, figsize=(15, 5))
        
        # 2D PSF
        im1 = axes[0].imshow(psf, cmap='viridis', extent=[-1, 1, -1, 1])
        axes[0].set_title('Double Gaussian PSF (2D)')
        axes[0].set_xlabel('X Position (μm)')
        axes[0].set_ylabel('Y Position (μm)')
        plt.colorbar(im1, ax=axes[0], label='PSF Value')
        
        # 1D cross-section
        center = psf.shape[0] // 2
        x_profile = X[center, :]
        y_profile = psf[center, :]
        axes[1].plot(x_profile, y_profile, 'b-', linewidth=2)
        axes[1].set_title('PSF Cross-Section')
        axes[1].set_xlabel('Position (μm)')
        axes[1].set_ylabel('PSF Value')
        axes[1].grid(True, alpha=0.3)
        
        # Radial profile
        r = np.sqrt(X**2 + Y**2)
        r_flat = r.flatten()
        psf_flat = psf.flatten()
        
        # Sort by radius
        sort_idx = np.argsort(r_flat)
        r_sorted = r_flat[sort_idx]
        psf_sorted = psf_flat[sort_idx]
        
        # Average over radial bins
        r_bins = np.linspace(0, np.max(r_sorted), 50)
        psf_radial = []
        for i in range(len(r_bins)-1):
            mask = (r_sorted >= r_bins[i]) & (r_sorted < r_bins[i+1])
            if np.any(mask):
                psf_radial.append(np.mean(psf_sorted[mask]))
            else:
                psf_radial.append(0)
        
        axes[2].plot(r_bins[:-1], psf_radial, 'r-', linewidth=2)
        axes[2].set_title('Radial Profile')
        axes[2].set_xlabel('Radius (μm)')
        axes[2].set_ylabel('PSF Value')
        axes[2].grid(True, alpha=0.3)
        
        plt.tight_layout()
        plt.show()

# Usage example
psf_calculator = DoubleGaussianPSF(grid_size=128, feature_size=2.0)
psf, X, Y = psf_calculator.calculate_psf(sigma1=0.03, sigma2=0.1, 
                                        amplitude1=0.8, amplitude2=0.2)
metrics = psf_calculator.calculate_metrics(psf)
psf_calculator.visualize_psf(psf, X, Y)

print(f"Peak intensity: {metrics['peak_intensity']:.3f}")
print(f"FWHM: {metrics['fwhm']:.3f} μm")
print(f"Total energy: {metrics['total_energy']:.3f}")

import numpy as np
import matplotlib.pyplot as plt

class PartialCoherenceAnalyzer:
    def __init__(self, wavelength=193):
        self.wavelength = wavelength
        
    def calculate_coherence_length(self, spectral_width, sigma):
        """Calculate coherence length"""
        coherence_length = (self.wavelength**2 / (2 * np.pi * spectral_width)) * (1 - sigma**2)
        return coherence_length
    
    def calculate_resolution(self, na):
        """Calculate resolution using Rayleigh criterion"""
        resolution = 0.61 * self.wavelength / na
        return resolution
    
    def calculate_depth_of_focus(self, na):
        """Calculate depth of focus"""
        dof = self.wavelength / (2 * na**2)
        return dof
    
    def calculate_contrast(self, pitch, na, sigma):
        """Calculate contrast for given pitch"""
        normalized_pitch = pitch / (self.wavelength / na)
        contrast = np.exp(-sigma**2) * np.sin(np.pi * normalized_pitch) / (np.pi * normalized_pitch)
        return contrast
    
    def generate_pupil_function(self, na, pupil_cutoff, illumination_type):
        """Generate pupil function for different illumination types"""
        grid_size = 100
        x = np.linspace(-1, 1, grid_size)
        y = np.linspace(-1, 1, grid_size)
        X, Y = np.meshgrid(x, y)
        R = np.sqrt(X**2 + Y**2)
        
        pupil = np.zeros_like(R)
        
        if illumination_type == 'circular':
            pupil[R <= pupil_cutoff] = 1.0
        elif illumination_type == 'annular':
            inner_radius = 0.3
            pupil[(R >= inner_radius) & (R <= pupil_cutoff)] = 1.0
        elif illumination_type == 'quadrupole':
            angle = np.arctan2(Y, X)
            pupil[(R <= pupil_cutoff) & (np.abs(np.sin(2 * angle)) > 0.7)] = 1.0
        elif illumination_type == 'dipole':
            pupil[(R <= pupil_cutoff) & (np.abs(X) > np.abs(Y)) & (np.abs(X) > 0.3)] = 1.0
        
        return pupil, X, Y
    
    def analyze_coherence_effects(self, na_range, sigma_range, pitch_range):
        """Analyze coherence effects over parameter ranges"""
        results = {}
        
        # Coherence length analysis
        coherence_lengths = []
        for sigma in sigma_range:
            cl = self.calculate_coherence_length(0.5, sigma)
            coherence_lengths.append(cl)
        
        results['coherence_lengths'] = coherence_lengths
        
        # Resolution analysis
        resolutions = []
        for na in na_range:
            res = self.calculate_resolution(na)
            resolutions.append(res)
        
        results['resolutions'] = resolutions
        
        # Contrast analysis
        contrast_data = []
        for pitch in pitch_range:
            contrast = self.calculate_contrast(pitch, 0.85, 0.3)
            contrast_data.append(contrast)
        
        results['contrast_data'] = contrast_data
        
        return results
    
    def visualize_analysis(self, results, na_range, sigma_range, pitch_range):
        """Visualize coherence analysis results"""
        fig, axes = plt.subplots(2, 2, figsize=(12, 10))
        
        # Coherence length vs sigma
        axes[0, 0].plot(sigma_range, results['coherence_lengths'], 'b-o')
        axes[0, 0].set_xlabel('Coherence Factor (σ)')
        axes[0, 0].set_ylabel('Coherence Length (μm)')
        axes[0, 0].set_title('Coherence Length vs σ')
        axes[0, 0].grid(True, alpha=0.3)
        
        # Resolution vs NA
        axes[0, 1].plot(na_range, results['resolutions'], 'r-o')
        axes[0, 1].set_xlabel('Numerical Aperture')
        axes[0, 1].set_ylabel('Resolution (nm)')
        axes[0, 1].set_title('Resolution vs NA')
        axes[0, 1].grid(True, alpha=0.3)
        
        # Contrast vs pitch
        axes[1, 0].plot(pitch_range, results['contrast_data'], 'g-o')
        axes[1, 0].set_xlabel('Pitch (nm)')
        axes[1, 0].set_ylabel('Contrast')
        axes[1, 0].set_title('Contrast vs Pitch')
        axes[1, 0].grid(True, alpha=0.3)
        
        # Pupil functions
        illumination_types = ['circular', 'annular', 'quadrupole', 'dipole']
        for i, illum_type in enumerate(illumination_types):
            pupil, X, Y = self.generate_pupil_function(0.85, 0.9, illum_type)
            im = axes[1, 1].imshow(pupil, cmap='viridis', extent=[-1, 1, -1, 1])
            axes[1, 1].set_title(f'Pupil Function - {illum_type.title()}')
            axes[1, 1].set_xlabel('Normalized X')
            axes[1, 1].set_ylabel('Normalized Y')
            break  # Show only one for space
        
        plt.tight_layout()
        plt.show()

# Usage example
analyzer = PartialCoherenceAnalyzer(wavelength=193)

# Parameter ranges
na_range = np.linspace(0.5, 1.0, 20)
sigma_range = np.linspace(0.1, 1.0, 20)
pitch_range = np.linspace(50, 500, 50)

# Run analysis
results = analyzer.analyze_coherence_effects(na_range, sigma_range, pitch_range)
analyzer.visualize_analysis(results, na_range, sigma_range, pitch_range)

# Calculate specific values
coherence_length = analyzer.calculate_coherence_length(0.5, 0.3)
resolution = analyzer.calculate_resolution(0.85)
dof = analyzer.calculate_depth_of_focus(0.85)

print(f"Coherence length: {coherence_length:.2f} μm")
print(f"Resolution: {resolution:.1f} nm")
print(f"Depth of focus: {dof:.2f} μm")

# Complete DUV Analysis Workflow
# This notebook demonstrates the complete DUV analysis workflow

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from duv_simulator import DUVSimulator, SimulationConfig

# 1. Initialize simulator
simulator = DUVSimulator()

# 2. Configure simulation parameters
config = SimulationConfig(
    wavelength=193,
    na=0.85,
    coherence_factor=0.3,
    resist_thickness=200.0,
    refractive_index=1.7
)

# 3. Run Monte Carlo simulation
print("Running Monte Carlo simulation...")
mc_results = simulator.run_monte_carlo(
    num_particles=100000,
    wavelength=193,
    na=0.85,
    feature_size=50.0,
    dose=30.0
)

# 4. Calculate Double Gaussian PSF
print("Calculating Double Gaussian PSF...")
psf_results = simulator.calculate_double_gaussian(
    sigma1=30.0,
    sigma2=100.0,
    amplitude1=0.8,
    amplitude2=0.2
)

# 5. Analyze partial coherence
print("Analyzing partial coherence...")
coherence_results = simulator.analyze_partial_coherence(
    na=0.85,
    sigma=0.3,
    wavelength=193
)

# 6. Model flare effects
print("Modeling flare effects...")
flare_results = simulator.model_flare(
    flare_level=0.02,
    flare_sigma=10.0,
    main_sigma=50.0
)

# 7. Calculate swing curves
print("Calculating swing curves...")
swing_results = simulator.calculate_swing_curves(
    wavelength=193,
    na=0.85,
    analysis_type="pitch"
)

# 8. Compare models
print("Comparing models...")
comparison_results = simulator.compare_models(
    feature_size=50.0,
    wavelength=193,
    na=0.85
)

# 9. Generate comprehensive report
print("Generating report...")
report = simulator.generate_report(
    mc_results=mc_results,
    psf_results=psf_results,
    coherence_results=coherence_results,
    flare_results=flare_results,
    swing_results=swing_results,
    comparison_results=comparison_results,
    include_plots=True,
    export_format="pdf"
)

# 10. Display results
print("\\n=== DUV Analysis Results ===")
print(f"Monte Carlo - Total Energy: {mc_results.statistics.total_energy:.2e} J")
print(f"Monte Carlo - Peak Intensity: {mc_results.statistics.peak_intensity:.2e} W/cm²")
print(f"PSF - FWHM1: {psf_results.metrics.fwhm1:.1f} nm")
print(f"PSF - FWHM2: {psf_results.metrics.fwhm2:.1f} nm")
print(f"Coherence - Coherence Length: {coherence_results.coherence_length:.2f} μm")
print(f"Coherence - Resolution: {coherence_results.resolution:.1f} nm")
print(f"Flare - Flare Intensity: {flare_results.metrics.flare_intensity:.1%}")
print(f"Flare - Contrast Reduction: {flare_results.metrics.contrast_reduction:.1%}")
print(f"Swing - Max Contrast: {swing_results.metrics.max_contrast:.3f}")
print(f"Swing - Max NILS: {swing_results.metrics.max_nils:.3f}")
print(f"Comparison - MC Accuracy: {comparison_results.monte_carlo.accuracy:.1f}%")
print(f"Comparison - DG Accuracy: {comparison_results.double_gaussian.accuracy:.1f}%")

# 11. Create summary plots
fig, axes = plt.subplots(2, 3, figsize=(18, 12))

# Monte Carlo energy deposition
axes[0, 0].imshow(mc_results.energy_deposition, cmap='viridis')
axes[0, 0].set_title('Monte Carlo Energy Deposition')
axes[0, 0].set_xlabel('X Position (nm)')
axes[0, 0].set_ylabel('Y Position (nm)')

# Double Gaussian PSF
axes[0, 1].imshow(psf_results.psf_2d, cmap='viridis')
axes[0, 1].set_title('Double Gaussian PSF')
axes[0, 1].set_xlabel('X Position (nm)')
axes[0, 1].set_ylabel('Y Position (nm)')

# Partial coherence pupil function
axes[0, 2].imshow(coherence_results.pupil_function, cmap='viridis')
axes[0, 2].set_title('Pupil Function')
axes[0, 2].set_xlabel('Normalized X')
axes[0, 2].set_ylabel('Normalized Y')

# Flare profile
axes[1, 0].plot(flare_results.profile['x'], flare_results.profile['main'], 'b-', label='Main PSF')
axes[1, 0].plot(flare_results.profile['x'], flare_results.profile['flare'], 'r-', label='Flare PSF')
axes[1, 0].plot(flare_results.profile['x'], flare_results.profile['total'], 'g--', label='Total PSF')
axes[1, 0].set_title('Flare Profile')
axes[1, 0].set_xlabel('Position (μm)')
axes[1, 0].set_ylabel('PSF Value')
axes[1, 0].legend()

# Swing curves
axes[1, 1].plot(swing_results.pitch_data, swing_results.contrast_data, 'b-o')
axes[1, 1].set_title('Swing Curves')
axes[1, 1].set_xlabel('Pitch (nm)')
axes[1, 1].set_ylabel('Contrast')
axes[1, 1].grid(True, alpha=0.3)

# Model comparison
comparison_metrics = ['Accuracy', 'Speed', 'Memory', 'Precision']
mc_scores = [comparison_results.monte_carlo.accuracy, 
            100 - comparison_results.monte_carlo.speed * 10,
            100 - comparison_results.monte_carlo.memory,
            comparison_results.monte_carlo.precision * 100]
dg_scores = [comparison_results.double_gaussian.accuracy,
            100 - comparison_results.double_gaussian.speed * 10,
            100 - comparison_results.double_gaussian.memory,
            comparison_results.double_gaussian.precision * 100]

x = np.arange(len(comparison_metrics))
width = 0.35

axes[1, 2].bar(x - width/2, mc_scores, width, label='Monte Carlo', color='blue', alpha=0.7)
axes[1, 2].bar(x + width/2, dg_scores, width, label='Double Gaussian', color='red', alpha=0.7)
axes[1, 2].set_title('Model Comparison')
axes[1, 2].set_xlabel('Metrics')
axes[1, 2].set_ylabel('Score')
axes[1, 2].set_xticks(x)
axes[1, 2].set_xticklabels(comparison_metrics)
axes[1, 2].legend()

plt.tight_layout()
plt.show()

print("\\nAnalysis complete! Report saved as 'duv_analysis_report.pdf'")