type ChannelBandType string

const (
    BandLongPeriodHigh    ChannelBandType = "L"  // 1 Hz
    BandLongPeriodLow     ChannelBandType = "L"  // 1 Hz  
    BandBroadband         ChannelBandType = "B"  // 10-100 Hz
    BandShortPeriodHigh  ChannelBandType = "H"  // 10-100 Hz
    BandShortPeriodLow   ChannelBandType = "S"  // 10-100 Hz
    BandMidPeriod        ChannelBandType = "M"  // 1-10 Hz
    BandVeryLongPeriod   ChannelBandType = "V"  // 0.01-0.1 Hz
)

type ChannelInstrumentType string

const (
    InstrumentHighGain       ChannelInstrumentType = "H"  // High gain seismometer
    InstrumentLowGain        ChannelInstrumentType = "L"  // Low gain seismometer
    InstrumentBroadband      ChannelInstrumentType = "B"  // Broadband seismometer
    InstrumentGeophone       ChannelInstrumentType = "G"  // Geophone
    InstrumentAccelerometer  ChannelInstrumentType = "N"  // Accelerometer
)

type ChannelOrientationType string

const (
    OrientationVertical   ChannelOrientationType = "Z"  // Vertical
    OrientationNorth      ChannelOrientationType = "N"  // North
    OrientationEast       ChannelOrientationType = "E"  // East
    Orientation1           ChannelOrientationType = "1"  // Orthogonal 1
    Orientation2           ChannelOrientationType = "2"  // Orthogonal 2
    Orientation3           ChannelOrientationType = "3"  // Orthogonal 3
)

// WebSocket handler for real-time streaming
func (h *Handler) handleWebSocket(c *gin.Context) {
    upgrader := websocket.Upgrader{
        CheckOrigin: func(r *http.Request) bool { return true },
    }
    
    conn, err := upgrader.Upgrade(c.Writer, c.Request, nil)
    if err != nil {
        return
    }
    defer conn.Close()
    
    // Stream waveform data
    for event := range h.events {
        conn.WriteJSON(event)
    }
}

// IASP91 velocity model (simplified)
type IASP91Model struct {
    // P-wave velocity layers (km/s)
    PLayers []VelocityLayer
    // S-wave velocity layers (km/s)
    SLayers []VelocityLayer
}

// Travel time calculation
func (m *IASP91Model) TravelTime(phase string, distance, depth float64) float64 {
    // Interpolate travel time from tables
    // Apply depth correction
    // Apply ellipticity correction
    return baseTime + depthCorrection + ellipticityCorrection
}

// STA/LTA ratio calculation
func (p *PhasePicker) Pick(samples []float64, sampleRate float64) []Pick {
    shortWindow := int(p.ShortTermWindow * sampleRate)
    longWindow := int(p.LongTermWindow * sampleRate)
    
    sta := calculateSTA(samples, shortWindow)
    lta := calculateLTA(samples, longWindow)
    
    ratio := sta / lta
    
    // Find triggers
    triggers := findTriggers(ratio, p.TriggerThreshold)
    
    // Refine with AIC
    return refineWithAIC(triggers, samples)
}

// Iterative Gauss-Newton location
func (g *Geiger) Locate(arrivals []PhaseArrival) Location {
    location := g.initialEstimate(arrivals)
    
    for i := 0; i < g.MaxIterations; i++ {
        // Calculate residuals
        residuals := g.calculateResiduals(arrivals, location)
        
        // Build Jacobian matrix
        J := g.buildJacobian(arrivals, location)
        
        // Solve normal equations: J^T * J * Δx = J^T * residuals
        delta := g.solveNormal(J, residuals)
        
        // Update location
        location.Lat += delta.Lat
        location.Lon += delta.Lon
        location.Depth += delta.Depth
        location.Time += delta.Time
        
        // Check convergence
        if g.converged(delta, g.ConvergenceThreshold) {
            break
        }
    }
    
    return location
}

func (l *ConfigLoader) LoadAll() error {
    return filepath.WalkDir(l.configDir, func(path string, d fs.DirEntry, err error) error {
        if d.IsDir() || !strings.HasSuffix(path, ".json") {
            return nil
        }
        return l.loadFile(path)
    })
}

// NLNM power in dB at frequency f
func (n *PetersonNLNM) Power(freq float64) float64 {
    // Peterson 1993 coefficients
    return interpolate(n.frequencies, n.powers, freq)
}

// Frequencies: 0.01, 0.1, 1, 10 Hz
// Powers: -187, -166, -166, -166 dB

func (g *PinkNoiseGenerator) Generate(n int) []float64 {
    samples := make([]float64, n)
    for i := 0; i < n; i++ {
        samples[i] = g.nextSample()
    }
    return samples
}

func (g *PinkNoiseGenerator) nextSample() float64 {
    // Voss-McCartney algorithm for 1/f noise
    k := g.key
    g.key++
    
    changed := k ^ g.key
    var sum float64
    for i := 0; i < 5; i++ {
        if changed&(1<<i) != 0 {
            g.white[i] = rand.NormFloat64()
        }
        sum += g.white[i]
    }
    sum += rand.NormFloat64()
    return sum / 6.0
}

type PAZResponse struct {
    Poles   []complex128  // Complex poles
    Zeros   []complex128  // Complex zeros
    Gain    float64       // Instrument gain
    A0      float64        // Normalization factor
}

// Calculate response at frequency
func (r *PAZResponse) Response(f float64) complex128 {
    s := complex(0, 2*math.Pi*f)
    
    // Numerator: zeros
    num := complex(1, 0)
    for _, z := range r.Zeros {
        num *= (s - z)
    }
    
    // Denominator: poles
    den := complex(1, 0)
    for _, p := range r.Poles {
        den *= (s - p)
    }
    
    return r.A0 * r.Gain * num / den
}

func TestEndToEndWorkflow(t *testing.T) {
    // 1. Station Generation
    station := generateStation("ASAR")
    
    // 2. Waveform Generation
    waveforms := generateWaveforms(station, distance, 40.0)
    
    // 3. Noise Model
    noise := generateNoise(PetersonNLNM)
    
    // 4. Channel Response
    response := calculateResponse(instrument, frequency)
    
    // 5. Travel Time
    travelTime := calculateTravelTime(distance, depth)
    
    // 6. Phase Picking
    picks := detectPhases(waveforms, 40.0)
    
    // 7. FK Analysis
    fk := analyzeFK(waveforms, array)
    
    // 8. Event Location
    location := locateEvent(picks, travelTimes)
    
    // 9. Magnitude
    magnitude := calculateMagnitude(amplitudes, distances)
    
    // 10. Validation
    result := validate(station, waveforms, picks, location, magnitude)
    
    // Assert accuracy
    if result.OverallAccuracy < 89.0 {
        t.Errorf("Accuracy %.1f%% below 89%%", result.OverallAccuracy)
    }
}

func (v *RealWaveformValidator) calculateNoiseMatch(simulated, expected []float64) float64 {
    // 1. Compare noise levels in pre-signal portion
    noiseWindow := min(len(simulated)/10, 100)
    
    simNoiseMean := mean(simulated[:noiseWindow])
    simNoiseStd := stdDev(simulated[:noiseWindow])
    expNoiseMean := mean(expected[:noiseWindow])
    expNoiseStd := stdDev(expected[:noiseWindow])
    
    // 2. Standard deviation ratio
    stdRatio := simNoiseStd / expNoiseStd
    if stdRatio > 1 {
        stdRatio = 1.0 / stdRatio
    }
    
    // 3. Mean comparison
    meanDiff := math.Abs(simNoiseMean - expNoiseMean)
    meanTolerance := (expNoiseStd + simNoiseStd) / 2
    meanAcc := 100.0 * (1.0 - meanDiff/meanTolerance)
    
    // 4. Power spectral density comparison
    psdMatch := comparePowerSpectralDensity(simulated, expected)
    
    // 5. Combined accuracy
    return 0.5 * (stdRatio * 100 + meanAcc) + 0.5 * psdMatch
}

func (v *RealWaveformValidator) calculateTimingAccuracy(simulated, expected []float64) float64 {
    // 1. Cross-correlation
    correlation := crossCorrelate(simulated, expected)
    correlationAcc := (correlation + 1.0) * 50.0
    
    // 2. Peak timing comparison
    simPeakIdx := findPeakIndex(simulated)
    expPeakIdx := findPeakIndex(expected)
    
    sampleTimingError := math.Abs(float64(simPeakIdx - expPeakIdx))
    
    timingFromPeaks := 100.0
    if sampleTimingError > 1 {
        timingFromPeaks = math.Max(0, 100.0 - (sampleTimingError-1)*5)
    }
    
    // 3. Combined: 50% correlation + 50% peak timing
    return 0.5 * correlationAcc + 0.5 * timingFromPeaks
}

func (c *Calibrator) CalibrateAmplitude(observed, expected float64, instrument string, freq float64) float64 {
    // Get instrument response
    cal := c.calibrations[instrument]
    
    // Calculate response at frequency
    response := c.frequencyResponse(cal, freq)
    
    // Correct amplitude
    return observed / response
}

// Full channel name: BHZ
// B = Broadband (10-100 Hz)
// H = High gain seismometer
// Z = Vertical component

func ParseChannelName(name string) (band, instrument, orientation string) {
    if len(name) >= 3 {
        band = string(name[0])
        instrument = string(name[1])
        orientation = string(name[2])
    }
    return
}

func (m *IASP91Model) interpolateTravelTime(phase string, distance, depth float64) float64 {
    // Find surrounding points in table
    d1, d2 := findSurrounding(distance)
    t1, t2 := findSurroundingTimes(phase, distance)
    
    // Bilinear interpolation
    t := t1 + (t2-t1)*(distance-d1)/(d2-d1)
    
    // Depth correction
    t += depth * m.depthCorrection(phase)
    
    return t
}

-- Stations
CREATE TABLE station (
    id VARCHAR(20) PRIMARY KEY,
    name VARCHAR(100),
    network VARCHAR(10),
    latitude DOUBLE PRECISION,
    longitude DOUBLE PRECISION,
    elevation DOUBLE PRECISION,
    description TEXT
);

-- Channels
CREATE TABLE channel (
    id VARCHAR(20) PRIMARY KEY,
    station_id VARCHAR(20) REFERENCES station(id),
    channel_type VARCHAR(10),  -- BHZ, BH1, BH2, etc.
    sample_rate DOUBLE PRECISION,
    depth DOUBLE PRECISION,
    azimuth DOUBLE PRECISION,
    dip DOUBLE PRECISION
);

-- Channel Groups
CREATE TABLE channel_group (
    id VARCHAR(20) PRIMARY KEY,
    name VARCHAR(100),
    station_id VARCHAR(20),
    channels TEXT[]  -- Array of channel IDs
);

-- FK Spectra Definitions
CREATE TABLE fk_spectra_definition (
    id VARCHAR(50) PRIMARY KEY,
    name VARCHAR(100),
    station_id VARCHAR(20),
    channel_group VARCHAR(50),
    phase VARCHAR(10),
    slow_start_x DOUBLE PRECISION,
    slow_start_y DOUBLE PRECISION,
    slow_delta_x DOUBLE PRECISION,
    slow_delta_y DOUBLE PRECISION,
    slow_count_x INTEGER,
    slow_count_y INTEGER,
    window_length DOUBLE PRECISION,
    sample_rate DOUBLE PRECISION
);

-- Beamforming Configuration
CREATE TABLE beamforming_configuration (
    id VARCHAR(50) PRIMARY KEY,
    name VARCHAR(100),
    station_id VARCHAR(20),
    channels TEXT[],
    beam_type VARCHAR(20),
    phase VARCHAR(10)
);

// Voss-McCartney pink noise
func GeneratePinkNoise(n int) []float64 {
    b := make([]float64, 5)  // 5 octaves
    result := make([]float64, n)
    
    for i := 0; i < n; i++ {
        // Update least significant changed bit
        k := i ^ (i - 1)
        
        for j := 0; j < 5; j++ {
            if k & (1 << j) != 0 {
                b[j] = rand.NormFloat64()
            }
        }
        
        // Sum octaves
        sum := 0.0
        for _, v := range b {
            sum += v
        }
        
        result[i] = sum / math.Sqrt(5.0)
    }
    
    return result
}

// Short-Term Average / Long-Term Average
func (p *PhasePicker) Detect(samples []float64, sampleRate float64) []PhasePick {
    shortWindow := int(p.ShortTerm * sampleRate)  // 1 second
    longWindow := int(p.LongTerm * sampleRate)     // 10 seconds
    
    sta := make([]float64, len(samples))
    lta := make([]float64, len(samples))
    
    // Calculate STA and LTA
    for i := longWindow; i < len(samples); i++ {
        // Short-term average
        staSum := 0.0
        for j := i - shortWindow; j < i; j++ {
            staSum += samples[j] * samples[j]
        }
        sta[i] = math.Sqrt(staSum / float64(shortWindow))
        
        // Long-term average
        ltaSum := 0.0
        for j := i - longWindow; j < i; j++ {
            ltaSum += samples[j] * samples[j]
        }
        lta[i] = math.Sqrt(ltaSum / float64(longWindow))
    }
    
    // Find triggers
    var picks []PhasePick
    triggered := false
    
    for i := longWindow; i < len(samples); i++ {
        if lta[i] > 0 {
            ratio := sta[i] / lta[i]
            
            if !triggered && ratio > p.Threshold {
                triggered = true
                triggerTime := float64(i) / sampleRate
                
                // Refine with AIC
                refinedTime := refineWithAIC(samples, i, sampleRate)
                
                picks = append(picks, PhasePick{
                    Time: refinedTime,
                    Phase: "P",
                    Amplitude: math.Abs(samples[i]),
                    SNR: ratio,
                    Confidence: calculateConfidence(ratio, p.Threshold),
                })
            }
            
            if triggered && ratio < p.Threshold * p.DetriggerRatio {
                triggered = false
            }
        }
    }
    
    return picks
}

// Frequency-Wavenumber Analysis
func AnalyzeFK(waveforms map[string][]float64, stations []Station, params FKParams) FKResult {
    n := len(stations)
    freqBins := params.FrequencyBins
    slowBins := params.SlownessBins
    
    // Compute cross-spectral matrix
    S := make([][]complex128, n)
    for i := 0; i < n; i++ {
        S[i] = make([]complex128, n)
        for j := 0; j < n; j++ {
            S[i][j] = crossSpectrum(waveforms[i], waveforms[j], freqBins)
        }
    }
    
    // Beamform over slowness grid
    maxPower := 0.0
    bestSlowX := 0.0
    bestSlowY := 0.0
    
    for sx := slowBins.MinX; sx <= slowBins.MaxX; sx += slowBins.DeltaX {
        for sy := slowBins.MinY; sy <= slowBins.MaxY; sy += slowBins.DeltaY {
            // Steering vector
            v := computeSteeringVector(stations, sx, sy, freqBins)
            
            // Beam power
            power := beamPower(S, v)
            
            if power > maxPower {
                maxPower = power
                bestSlowX = sx
                bestSlowY = sy
            }
        }
    }
    
    return FKResult{
        Power: maxPower,
        SlownessX: bestSlowX,
        SlownessY: bestSlowY,
        Azimuth: math.Atan2(bestSlowY, bestSlowX) * 180 / math.Pi,
    }
}

func (g *Geiger) Locate(arrivals []PhaseArrival) Location {
    // Initial estimate
    loc := g.initialEstimate(arrivals)
    
    for iter := 0; iter < g.MaxIterations; iter++ {
        // Calculate residuals
        residuals := make([]float64, len(arrivals))
        for i, arr := range arrivals {
            predictedTime := g.calculateTravelTime(arr.Station, loc)
            residuals[i] = arr.Time - predictedTime
        }
        
        // Build Jacobian matrix
        J := g.buildJacobian(arrivals, loc)
        
        // Solve: J^T * W * J * Δ = J^T * W * r
        delta := g.solveNormal(J, residuals)
        
        // Update location
        loc.Latitude += delta[0]
        loc.Longitude += delta[1]
        loc.Depth += delta[2]
        loc.OriginTime += delta[3]
        
        // Check convergence
        if g.converged(delta, g.Threshold) {
            break
        }
    }
    
    return loc
}

# k8s.yaml (abbreviated)
apiVersion: apps/v1
kind: Deployment
metadata:
  name: dod-simulator
  namespace: gms
spec:
  replicas: 1
  selector:
    matchLabels:
      app: dod-simulator
  template:
    spec:
      containers:
      - name: dod-simulator
        image: dod-simulator:latest
        command: ["/app/dod-simulator", "serve", "--host", "0.0.0.0", "--port", "8081"]
        ports:
        - containerPort: 8081
        env:
        - name: DOD_AUTH_ENABLED
          value: "true"

// Complete SEED band types
type ChannelBandType string

const (
    BandLongPeriodHigh   ChannelBandType = "L"  // 1 Hz, High gain
    BandLongPeriodLow    ChannelBandType = "L"  // 1 Hz, Low gain
    BandBroadband        ChannelBandType = "B"  // 10-100 Hz
    BandShortPeriodHigh  ChannelBandType = "H"  // 10-100 Hz, High gain
    BandShortPeriodLow   ChannelBandType = "S"  // 10-100 Hz, Low gain
    BandMidPeriod        ChannelBandType = "M"  // 1-10 Hz
    BandVeryLongPeriod   ChannelBandType = "V"  // 0.01-0.1 Hz
)

// Instrument types
type ChannelInstrumentType string

const (
    InstrumentHighGain      ChannelInstrumentType = "H"
    InstrumentLowGain       ChannelInstrumentType = "L"
    InstrumentBroadband     ChannelInstrumentType = "B"
    InstrumentGeophone      ChannelInstrumentType = "G"
    InstrumentAccelerometer ChannelInstrumentType = "N"
    InstrumentGeophoneHigh  ChannelInstrumentType = "D"
    InstrumentGeophoneLow  ChannelInstrumentType = "E"
    InstrumentHydrophone    ChannelInstrumentType = "J"
    InstrumentDipole        ChannelInstrumentType = "M"
)

// Orientation types
type ChannelOrientationType string

const (
    OrientationVertical  ChannelOrientationType = "Z"
    OrientationNorth     ChannelOrientationType = "N"
    OrientationEast      ChannelOrientationType = "E"
    Orientation1          ChannelOrientationType = "1"
    Orientation2         ChannelOrientationType = "2"
    Orientation3         ChannelOrientationType = "3"
)

// IASP91 velocity model
type IASP91Model struct {
    // P-wave velocity profile (km/s)
    PLayers []VelocityLayer
    // S-wave velocity profile (km/s)
    SLayers []VelocityLayer
}

type VelocityLayer struct {
    Depth      float64  // km
    Velocity   float64  // km/s
}

// Travel time calculation
func (m *IASP91Model) TravelTime(phase string, distance, depth float64) float64 {
    // Interpolate from IASP91 tables
    // Apply depth correction
    // Apply ellipticity correction
    // Return travel time in seconds
}

// Ellipticity correction
func (m *IASP91Model) EllipticityCorrection(stationLat, eventLat, time, phase string) float64 {
    // Based on Dziewonski & Gilbert (1976)
    // Up to 0.5 seconds correction
}

// Basin correction
func (m *IASP91Model) BasinCorrection(stationLon, stationLat float64) float64 {
    // Major basin corrections
    basins := map[string]float64{
        "LA": 0.15,    // Los Angeles Basin
        "SF": 0.12,    // San Francisco Bay
        "NY": 0.08,    // New York Basin
    }
    // Return basin-specific correction
}

// Peterson New Low Noise Model (NLNM)
type PetersonNLNM struct {
    frequencies []float64
    powers      []float64
}

func NewPetersonNLNM() *PetersonNLNM {
    return &PetersonNLNM{
        frequencies: []float64{
            0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 1, 2, 5, 10,
        },
        powers: []float64{
            -187.0, -182.0, -175.0, -166.0, -159.0, -151.0, -144.0, -139.0, -127.0, -113.0,
        },
    }
}

func (n *PetersonNLNM) GetPowerDB(freq float64) float64 {
    // Log-linear interpolation
    return interpolateLog(n.frequencies, n.powers, freq)
}

// Peterson New High Noise Model (NHNM)
type PetersonNHNM struct {
    frequencies []float64
    powers      []float64
}

func NewPetersonNHNM() *PetersonNHNM {
    return &PetersonNHNM{
        frequencies: []float64{
            0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 1, 2, 5, 10,
        },
        powers: []float64{
            -108.0, -108.0, -108.0, -108.0, -108.0, -105.0, -102.0, -97.0, -87.0, -73.0,
        },
    }
}

// Phase-specific parameters
type PhaseSpecificParams struct {
    ShortTermWindow   float64
    LongTermWindow    float64
    TriggerThreshold  float64
    DetriggerRatio    float64
    MinPeakAmplitude  float64
    MinSNR            float64
}

var (
    PWaveParams = PhaseSpecificParams{
        ShortTermWindow:   1.0,
        LongTermWindow:    10.0,
        TriggerThreshold:  3.0,
        DetriggerRatio:    0.5,
        MinPeakAmplitude:  0.5,
        MinSNR:            2.0,
    }
    
    SWaveParams = PhaseSpecificParams{
        ShortTermWindow:   2.0,
        LongTermWindow:    15.0,
        TriggerThreshold:  2.5,
        DetriggerRatio:    0.6,
        MinPeakAmplitude:  0.4,
        MinSNR:            1.8,
    }
    
    SurfaceWaveParams = PhaseSpecificParams{
        ShortTermWindow:   5.0,
        LongTermWindow:    30.0,
        TriggerThreshold:  2.0,
        DetriggerRatio:    0.7,
        MinPeakAmplitude:  0.3,
        MinSNR:            1.5,
    }
)

// Geiger iterative location
func (g *Geiger) Locate(arrivals []PhaseArrival) Location {
    // Initial estimate from first arrivals
    loc := g.initialEstimate(arrivals)
    
    damping := g.Damping
    
    for iter := 0; iter < g.MaxIterations; iter++ {
        // Calculate residuals
        residuals := make([]float64, len(arrivals))
        for i, arr := range arrivals {
            predictedTime := g.calculateTravelTime(arr.Station, loc)
            residuals[i] = arr.Time - predictedTime
        }
        
        // Build Jacobian matrix
        J := make([][]float64, len(arrivals))
        for i := range J {
            J[i] = make([]float64, 4)  // lat, lon, depth, time
            J[i][0] = g.dTravelTime_dLat(arr.Station, loc)
            J[i][1] = g.dTravelTime_dLon(arr.Station, loc)
            J[i][2] = g.dTravelTime_dDepth(arr.Station, loc)
            J[i][3] = 1.0  // time derivative = 1
        }
        
        // Apply damping
        for i := 0; i < 4; i++ {
            J[i][i] += damping
        }
        
        // Solve normal equations
        delta := g.solveNormal(J, residuals)
        
        // Update location
        loc.Latitude += delta[0]
        loc.Longitude += delta[1]
        loc.Depth += delta[2]
        loc.OriginTime += delta[3]
        
        // Check convergence
        if g.converged(delta, g.Threshold) {
            break
        }
    }
    
    return loc
}

// FK analysis for beamforming
func AnalyzeFK(waveforms map[string][]float64, stations []Station, params FKParams) FKResult {
    n := len(stations)
    frequencies := params.Frequencies
    
    // Compute cross-spectral matrix
    S := make([][]complex128, n)
    for i := 0; i < n; i++ {
        S[i] = make([]complex128, n)
        for j := 0; j < n; j++ {
            S[i][j] = crossSpectrum(waveforms[i], waveforms[j], frequencies)
        }
    }
    
    // Search slowness space
    var maxPower float64
    var bestSx, bestSy float64
    
    for sx := params.SlowMinX; sx <= params.SlowMaxX; sx += params.SlowDeltaX {
        for sy := params.SlowMinY; sy <= params.SlowMaxY; sy += params.SlowDeltaY {
            // Compute beam power
            power := beamPower(S, stations, sx, sy, frequencies)
            
            if power > maxPower {
                maxPower = power
                bestSx = sx
                bestSy = sy
            }
        }
    }
    
    // Compute azimuth from slowness
    azimuth := math.Atan2(bestSy, bestSx) * 180 / math.Pi
    if azimuth < 0 {
        azimuth += 360
    }
    
    return FKResult{
        Power:      maxPower,
        SlownessX:  bestSx,
        SlownessY:  bestSy,
        Slowness:   math.Sqrt(bestSx*bestSx + bestSy*bestSy),
        Azimuth:    azimuth,
    }
}

// WebSocket streaming handler
func (h *Handler) handleWebSocket(c *gin.Context) {
    upgrader := websocket.Upgrader{
        CheckOrigin: func(r *http.Request) bool {
            return true // Allow all origins
        },
    }
    
    conn, err := upgrader.Upgrade(c.Writer, c.Request, nil)
    if err != nil {
        return
    }
    defer conn.Close()
    
    // Subscribe to channels
    channels := parseChannels(c.Query("channels"))
    h.streamManager.Subscribe(conn, channels)
    defer h.streamManager.Unsubscribe(conn)
    
    // Stream events
    for event := range h.eventChannel {
        data, _ := json.Marshal(event)
        conn.WriteMessage(websocket.TextMessage, data)
    }
}

// Stream event phases
func StreamEvent(sm *StreamManager, stationID string, channels []string, eventTime time.Time) {
    for _, channel := range channels {
        go func(ch string) {
            // Stream P-wave
            StreamEventPhase(sm, stationID, ch, "P", eventTime, 5*time.Second)
            // Stream S-wave
            StreamEventPhase(sm, stationID, ch, "S", eventTime.Add(5*time.Second), 8*time.Second)
            // Stream surface wave
            StreamEventPhase(sm, stationID, ch, "L", eventTime.Add(13*time.Second), 15*time.Second)
        }(channel)
    }
}

// Real waveform validator
type RealWaveformValidator struct {
    tolerance float64
}

func (v *RealWaveformValidator) Validate(simulated, expected []float64, sampleRate float64) ValidationResult {
    // Phase accuracy
    phaseAcc := v.calculatePhaseAccuracy(simulated, expected, sampleRate)
    
    // Amplitude accuracy
    ampAcc := v.calculateAmplitudeAccuracy(simulated, expected)
    
    // Timing accuracy
    timeAcc := v.calculateTimingAccuracy(simulated, expected)
    
    // Noise match
    noiseAcc := v.calculateNoiseMatch(simulated, expected)
    
    // Overall
    overall := (phaseAcc + ampAcc + timeAcc + noiseAcc) / 4.0
    
    return ValidationResult{
        PhaseAccuracy:     phaseAcc,
        AmplitudeAccuracy: ampAcc,
        TimingAccuracy:    timeAcc,
        NoiseMatch:       noiseAcc,
        OverallAccuracy:   overall,
    }
}

{
  "id": "ASAR",
  "name": "Alice Springs Array",
  "network": "AU",
  "latitude": -23.664,
  "longitude": 133.905,
  "elevation": 540.0,
  "channels": [
    {"id": "AS31", "type": "BHZ", "sample_rate": 40},
    {"id": "AS31", "type": "BH1", "sample_rate": 40},
    {"id": "AS31", "type": "BH2", "sample_rate": 40}
  ]
}

{
  "name": "ASAR_FK",
  "stationId": "ASAR",
  "channelGroupId": "AS31",
  "phase": "P",
  "slowStartX": -0.5,
  "slowStartY": -0.5,
  "slowDeltaX": 0.02,
  "slowDeltaY": 0.02,
  "slowCountX": 50,
  "slowCountY": 50,
  "windowLength": 10.0,
  "sampleRate": 40.0
}