Program Listing for File NuMagSANSlib_SpectralData.h

↰ Return to documentation for file (src/OutputData/SpectralData/NuMagSANSlib_SpectralData.h)

// File         : NuMagSANSlib_SANSData.h
// Author       : Dr. Michael Philipp ADAMS
// Company      : University of Luxembourg
// Department   : Department of Physics and Materials Sciences
// Group        : NanoMagnetism Group
// Group Leader : Prof. Andreas Michels
// Version      : 28 October 2025
// OS           : Linux Ubuntu
// Language     : CUDA C++

#include <iostream>
#include <fstream>
#include <sstream>
#include <sys/stat.h>
#include <sys/types.h>
#include <math.h>
#include <string>
#include <vector>
#include <stdlib.h>
#include <time.h>
#include <cuda_runtime.h>
#include <cublas_v2.h>
#include <stdexcept>
#include <math.h>
#include <chrono>
#include <dirent.h>
#include <unistd.h>
#include <filesystem>

#include "NuMagSANSlib_SpectralDataStruct.h"
#include "NuMagSANSlib_SpectralData2csv.h"
#include "NuMagSANSlib_SpectralData2hdf5.h"

using namespace std;

void scale_SpectralData(ScalingFactors* ScalFactors, SpectralData* SpecData, InputFileData* InputData) {

    unsigned int Nq = *SpecData->Nq;
    unsigned int k_max = *SpecData->k_max;
    unsigned int L = 2 * Nq * (k_max + 1);
    unsigned int Lk = 2 * (k_max + 1);

    for (unsigned long int l = 0; l < L; l++) {

        // basic SANS cross sections
        SpecData->I_Nuc_unpolarized[l] = SpecData->I_Nuc_unpolarized[l] * ScalFactors->Nuc_SANS_SF;
        SpecData->I_Mag_unpolarized[l] = SpecData->I_Mag_unpolarized[l] * ScalFactors->Mag_SANS_SF;
        SpecData->I_Mag_polarized[l] = SpecData->I_Mag_polarized[l] * ScalFactors->Mag_SANS_SF;
        ;
        SpecData->I_NucMag[l] = SpecData->I_NucMag[l] * ScalFactors->NucMag_SANS_SF;
        SpecData->I_Mag_chiral[l] = SpecData->I_Mag_chiral[l] * ScalFactors->Mag_SANS_SF;

        // composed SANS cross sections
        SpecData->I_Mag_spin_flip[l] = SpecData->I_Mag_unpolarized[l] - SpecData->I_Mag_polarized[l];
        SpecData->I_Mag_spin_flip_pm[l] = SpecData->I_Mag_spin_flip[l] + SpecData->I_Mag_chiral[l];
        SpecData->I_Mag_spin_flip_mp[l] = SpecData->I_Mag_spin_flip[l] - SpecData->I_Mag_chiral[l];
        SpecData->I_Mag_non_spin_flip_pp[l] =
            SpecData->I_Nuc_unpolarized[l] + SpecData->I_NucMag[l] + SpecData->I_Mag_polarized[l];
        SpecData->I_Mag_non_spin_flip_mm[l] =
            SpecData->I_Nuc_unpolarized[l] - SpecData->I_NucMag[l] + SpecData->I_Mag_polarized[l];
        SpecData->I_Mag_sanspol_p[l] = SpecData->I_Mag_non_spin_flip_pp[l] + SpecData->I_Mag_spin_flip_pm[l];
        SpecData->I_Mag_sanspol_m[l] = SpecData->I_Mag_non_spin_flip_mm[l] + SpecData->I_Mag_spin_flip_mp[l];
    }

    // =====================================================
    // (2) scaling and normalization of the amplitude spectrum
    // =====================================================
    auto normalizeAmplitude = [&](float* A_array, float scaleFactor) {
        // (a) Scaling
        for (unsigned int k = 0; k < Lk; ++k)
            A_array[k] *= scaleFactor;

        // (b) compute norm of the ampliutde spectrum
        float normSum = 0.0f;
        for (unsigned int k = 0; k <= k_max; ++k) {
            float Ac = A_array[k];
            float As = A_array[k + (k_max + 1)];
            normSum += sqrtf(Ac * Ac + As * As);
        }

        if (normSum == 0.0f)
            return;

        // (c) normalization
        for (unsigned int k = 0; k <= k_max; ++k) {
            A_array[k] /= normSum;
            A_array[k + (k_max + 1)] /= normSum;
        }
    };

    // normalizing of the base amplitudes
    normalizeAmplitude(SpecData->A_Nuc_unpolarized, ScalFactors->Nuc_SANS_SF);
    normalizeAmplitude(SpecData->A_Mag_unpolarized, ScalFactors->Mag_SANS_SF);
    normalizeAmplitude(SpecData->A_Mag_polarized, ScalFactors->Mag_SANS_SF);
    normalizeAmplitude(SpecData->A_NucMag, ScalFactors->NucMag_SANS_SF);
    normalizeAmplitude(SpecData->A_Mag_chiral, ScalFactors->Mag_SANS_SF);

    // =====================================================
    // composition of the secondary amplitudes
    // =====================================================
    for (unsigned int k = 0; k < Lk; ++k) {
        SpecData->A_Mag_spin_flip[k] = SpecData->A_Mag_unpolarized[k] - SpecData->A_Mag_polarized[k];
        SpecData->A_Mag_spin_flip_pm[k] = SpecData->A_Mag_spin_flip[k] + SpecData->A_Mag_chiral[k];
        SpecData->A_Mag_spin_flip_mp[k] = SpecData->A_Mag_spin_flip[k] - SpecData->A_Mag_chiral[k];
        SpecData->A_Mag_non_spin_flip_pp[k] =
            SpecData->A_Nuc_unpolarized[k] + SpecData->A_NucMag[k] + SpecData->A_Mag_polarized[k];
        SpecData->A_Mag_non_spin_flip_mm[k] =
            SpecData->A_Nuc_unpolarized[k] - SpecData->A_NucMag[k] + SpecData->A_Mag_polarized[k];
        SpecData->A_Mag_sanspol_p[k] = SpecData->A_Mag_non_spin_flip_pp[k] + SpecData->A_Mag_spin_flip_pm[k];
        SpecData->A_Mag_sanspol_m[k] = SpecData->A_Mag_non_spin_flip_mm[k] + SpecData->A_Mag_spin_flip_mp[k];
    }

    // normalize the secondary amplitudes
    normalizeAmplitude(SpecData->A_Mag_spin_flip, 1.0f);
    normalizeAmplitude(SpecData->A_Mag_spin_flip_pm, 1.0f);
    normalizeAmplitude(SpecData->A_Mag_spin_flip_mp, 1.0f);
    normalizeAmplitude(SpecData->A_Mag_non_spin_flip_pp, 1.0f);
    normalizeAmplitude(SpecData->A_Mag_non_spin_flip_mm, 1.0f);
    normalizeAmplitude(SpecData->A_Mag_sanspol_p, 1.0f);
    normalizeAmplitude(SpecData->A_Mag_sanspol_m, 1.0f);
}

void allocate_SpectralData_RAM(InputFileData* InputData, ScatteringData* SANSData, SpectralData* SpecData) {

    SpecData->Nq = (unsigned int*)malloc(sizeof(unsigned int));
    SpecData->Ntheta = (unsigned int*)malloc(sizeof(unsigned int));
    SpecData->k_max = (unsigned int*)malloc(sizeof(unsigned int));
    SpecData->dtheta = (float*)malloc(sizeof(float));

    *SpecData->Nq = InputData->N_q;
    *SpecData->Ntheta = InputData->N_theta;
    *SpecData->k_max = InputData->k_max;
    *SpecData->dtheta = *SANSData->dtheta;

    // allocate the q-space vector
    SpecData->q = (float*)malloc((*(SpecData->Nq)) * sizeof(float));
    memset(SpecData->q, 0, (*SpecData->Nq) * sizeof(float));
    memcpy(SpecData->q, SANSData->q_1D, (*SpecData->Nq) * sizeof(float));

    // allocate spectral arrays
    size_t len = 2 * (*(SpecData->Nq)) * ((*(SpecData->k_max)) + 1) * sizeof(float);
    SpecData->I_Nuc_unpolarized = (float*)malloc(len);
    SpecData->I_Mag_unpolarized = (float*)malloc(len);
    SpecData->I_Mag_polarized = (float*)malloc(len);
    SpecData->I_NucMag = (float*)malloc(len);
    SpecData->I_Mag_spin_flip = (float*)malloc(len);
    SpecData->I_Mag_chiral = (float*)malloc(len);
    SpecData->I_Mag_spin_flip_pm = (float*)malloc(len);
    SpecData->I_Mag_spin_flip_mp = (float*)malloc(len);
    SpecData->I_Mag_non_spin_flip_pp = (float*)malloc(len);
    SpecData->I_Mag_non_spin_flip_mm = (float*)malloc(len);
    SpecData->I_Mag_sanspol_p = (float*)malloc(len);
    SpecData->I_Mag_sanspol_m = (float*)malloc(len);

    // --- Initialize to zero ---
    memset(SpecData->I_Nuc_unpolarized, 0, len);
    memset(SpecData->I_Mag_unpolarized, 0, len);
    memset(SpecData->I_Mag_polarized, 0, len);
    memset(SpecData->I_NucMag, 0, len);
    memset(SpecData->I_Mag_spin_flip, 0, len);
    memset(SpecData->I_Mag_chiral, 0, len);
    memset(SpecData->I_Mag_spin_flip_pm, 0, len);
    memset(SpecData->I_Mag_spin_flip_mp, 0, len);
    memset(SpecData->I_Mag_non_spin_flip_pp, 0, len);
    memset(SpecData->I_Mag_non_spin_flip_mm, 0, len);
    memset(SpecData->I_Mag_sanspol_p, 0, len);
    memset(SpecData->I_Mag_sanspol_m, 0, len);

    // allocate spectral arrays
    size_t lenA = 2 * ((*(SpecData->k_max)) + 1) * sizeof(float);
    SpecData->A_Nuc_unpolarized = (float*)malloc(lenA);
    SpecData->A_Mag_unpolarized = (float*)malloc(lenA);
    SpecData->A_Mag_polarized = (float*)malloc(lenA);
    SpecData->A_NucMag = (float*)malloc(lenA);
    SpecData->A_Mag_spin_flip = (float*)malloc(lenA);
    SpecData->A_Mag_chiral = (float*)malloc(lenA);
    SpecData->A_Mag_spin_flip_pm = (float*)malloc(lenA);
    SpecData->A_Mag_spin_flip_mp = (float*)malloc(lenA);
    SpecData->A_Mag_non_spin_flip_pp = (float*)malloc(lenA);
    SpecData->A_Mag_non_spin_flip_mm = (float*)malloc(lenA);
    SpecData->A_Mag_sanspol_p = (float*)malloc(lenA);
    SpecData->A_Mag_sanspol_m = (float*)malloc(lenA);

    // --- Initialize to zero ---
    memset(SpecData->A_Nuc_unpolarized, 0, lenA);
    memset(SpecData->A_Mag_unpolarized, 0, lenA);
    memset(SpecData->A_Mag_polarized, 0, lenA);
    memset(SpecData->A_NucMag, 0, lenA);
    memset(SpecData->A_Mag_spin_flip, 0, lenA);
    memset(SpecData->A_Mag_chiral, 0, lenA);
    memset(SpecData->A_Mag_spin_flip_pm, 0, lenA);
    memset(SpecData->A_Mag_spin_flip_mp, 0, lenA);
    memset(SpecData->A_Mag_non_spin_flip_pp, 0, lenA);
    memset(SpecData->A_Mag_non_spin_flip_mm, 0, lenA);
    memset(SpecData->A_Mag_sanspol_p, 0, lenA);
    memset(SpecData->A_Mag_sanspol_m, 0, lenA);
}

void allocate_SpectralData_GPU(SpectralData* SpecData, SpectralData* SpecData_gpu) {

    // Allocate scalar members on GPU
    cudaMalloc(&(SpecData_gpu->Nq), sizeof(unsigned int));
    cudaMalloc(&(SpecData_gpu->Ntheta), sizeof(unsigned int));
    cudaMalloc(&(SpecData_gpu->k_max), sizeof(unsigned int));
    cudaMalloc(&(SpecData_gpu->dtheta), sizeof(float));

    // Copy scalar values from host SpecData → device SpecData_gpu
    cudaMemcpy(SpecData_gpu->Nq, SpecData->Nq, sizeof(unsigned int), cudaMemcpyHostToDevice);
    cudaMemcpy(SpecData_gpu->Ntheta, SpecData->Ntheta, sizeof(unsigned int), cudaMemcpyHostToDevice);
    cudaMemcpy(SpecData_gpu->k_max, SpecData->k_max, sizeof(unsigned int), cudaMemcpyHostToDevice);
    cudaMemcpy(SpecData_gpu->dtheta, SpecData->dtheta, sizeof(float), cudaMemcpyHostToDevice);

    unsigned int Nq = *SpecData->Nq;
    unsigned int k_max = *SpecData->k_max;
    size_t len_q = Nq * sizeof(float);
    size_t len_I = 2 * Nq * (k_max + 1) * sizeof(float);
    size_t len_A = 2 * (k_max + 1) * sizeof(float);

    // Allocate arrays on GPU
    cudaMalloc(&(SpecData_gpu->q), len_q);
    cudaMemcpy(SpecData_gpu->q, SpecData->q, len_q, cudaMemcpyHostToDevice);

    cudaMalloc(&(SpecData_gpu->I_Nuc_unpolarized), len_I);
    cudaMalloc(&(SpecData_gpu->I_Mag_unpolarized), len_I);
    cudaMalloc(&(SpecData_gpu->I_Mag_polarized), len_I);
    cudaMalloc(&(SpecData_gpu->I_NucMag), len_I);
    cudaMalloc(&(SpecData_gpu->I_Mag_spin_flip), len_I);
    cudaMalloc(&(SpecData_gpu->I_Mag_chiral), len_I);
    cudaMalloc(&(SpecData_gpu->I_Mag_spin_flip_pm), len_I);
    cudaMalloc(&(SpecData_gpu->I_Mag_spin_flip_mp), len_I);
    cudaMalloc(&(SpecData_gpu->I_Mag_non_spin_flip_pp), len_I);
    cudaMalloc(&(SpecData_gpu->I_Mag_non_spin_flip_mm), len_I);
    cudaMalloc(&(SpecData_gpu->I_Mag_sanspol_p), len_I);
    cudaMalloc(&(SpecData_gpu->I_Mag_sanspol_m), len_I);

    // Initialize GPU arrays with zeros
    cudaMemset(SpecData_gpu->I_Nuc_unpolarized, 0, len_I);
    cudaMemset(SpecData_gpu->I_Mag_unpolarized, 0, len_I);
    cudaMemset(SpecData_gpu->I_Mag_polarized, 0, len_I);
    cudaMemset(SpecData_gpu->I_NucMag, 0, len_I);
    cudaMemset(SpecData_gpu->I_Mag_spin_flip, 0, len_I);
    cudaMemset(SpecData_gpu->I_Mag_chiral, 0, len_I);
    cudaMemset(SpecData_gpu->I_Mag_spin_flip_pm, 0, len_I);
    cudaMemset(SpecData_gpu->I_Mag_spin_flip_mp, 0, len_I);
    cudaMemset(SpecData_gpu->I_Mag_non_spin_flip_pp, 0, len_I);
    cudaMemset(SpecData_gpu->I_Mag_non_spin_flip_mm, 0, len_I);
    cudaMemset(SpecData_gpu->I_Mag_sanspol_p, 0, len_I);
    cudaMemset(SpecData_gpu->I_Mag_sanspol_m, 0, len_I);

    cudaMalloc(&(SpecData_gpu->A_Nuc_unpolarized), len_A);
    cudaMalloc(&(SpecData_gpu->A_Mag_unpolarized), len_A);
    cudaMalloc(&(SpecData_gpu->A_Mag_polarized), len_A);
    cudaMalloc(&(SpecData_gpu->A_NucMag), len_A);
    cudaMalloc(&(SpecData_gpu->A_Mag_spin_flip), len_A);
    cudaMalloc(&(SpecData_gpu->A_Mag_chiral), len_A);
    cudaMalloc(&(SpecData_gpu->A_Mag_spin_flip_pm), len_A);
    cudaMalloc(&(SpecData_gpu->A_Mag_spin_flip_mp), len_A);
    cudaMalloc(&(SpecData_gpu->A_Mag_non_spin_flip_pp), len_A);
    cudaMalloc(&(SpecData_gpu->A_Mag_non_spin_flip_mm), len_A);
    cudaMalloc(&(SpecData_gpu->A_Mag_sanspol_p), len_A);
    cudaMalloc(&(SpecData_gpu->A_Mag_sanspol_m), len_A);

    // Initialize GPU arrays with zeros
    cudaMemset(SpecData_gpu->A_Nuc_unpolarized, 0, len_A);
    cudaMemset(SpecData_gpu->A_Mag_unpolarized, 0, len_A);
    cudaMemset(SpecData_gpu->A_Mag_polarized, 0, len_A);
    cudaMemset(SpecData_gpu->A_NucMag, 0, len_A);
    cudaMemset(SpecData_gpu->A_Mag_spin_flip, 0, len_A);
    cudaMemset(SpecData_gpu->A_Mag_chiral, 0, len_A);
    cudaMemset(SpecData_gpu->A_Mag_spin_flip_pm, 0, len_A);
    cudaMemset(SpecData_gpu->A_Mag_spin_flip_mp, 0, len_A);
    cudaMemset(SpecData_gpu->A_Mag_non_spin_flip_pp, 0, len_A);
    cudaMemset(SpecData_gpu->A_Mag_non_spin_flip_mm, 0, len_A);
    cudaMemset(SpecData_gpu->A_Mag_sanspol_p, 0, len_A);
    cudaMemset(SpecData_gpu->A_Mag_sanspol_m, 0, len_A);

    cudaMalloc(&SpecData_gpu, sizeof(SpectralData));
    cudaMemcpy(SpecData_gpu, SpecData, sizeof(SpectralData), cudaMemcpyHostToDevice);

    cudaDeviceSynchronize();
}

void init_SpectralData(InputFileData* InputData, ScatteringData* SANSData, SpectralData* SpecData,
                       SpectralData* SpecData_gpu) {

    allocate_SpectralData_RAM(InputData, SANSData, SpecData);
    allocate_SpectralData_GPU(SpecData, SpecData_gpu);
}

void copyGPU2RAM_SpectralData(SpectralData* S, SpectralData* S_gpu) {

    cudaDeviceSynchronize();

    unsigned int Nq = *S->Nq;
    unsigned int k_max = *S->k_max;

    size_t len_I = 2 * Nq * (k_max + 1) * sizeof(float);
    size_t len_A = 2 * (k_max + 1) * sizeof(float);

    // Copy spectral arrays
    cudaMemcpy(S->I_Nuc_unpolarized, S_gpu->I_Nuc_unpolarized, len_I, cudaMemcpyDeviceToHost);
    cudaMemcpy(S->I_Mag_unpolarized, S_gpu->I_Mag_unpolarized, len_I, cudaMemcpyDeviceToHost);
    cudaMemcpy(S->I_Mag_polarized, S_gpu->I_Mag_polarized, len_I, cudaMemcpyDeviceToHost);
    cudaMemcpy(S->I_NucMag, S_gpu->I_NucMag, len_I, cudaMemcpyDeviceToHost);
    cudaMemcpy(S->I_Mag_spin_flip, S_gpu->I_Mag_spin_flip, len_I, cudaMemcpyDeviceToHost);
    cudaMemcpy(S->I_Mag_chiral, S_gpu->I_Mag_chiral, len_I, cudaMemcpyDeviceToHost);
    cudaMemcpy(S->I_Mag_spin_flip_pm, S_gpu->I_Mag_spin_flip_pm, len_I, cudaMemcpyDeviceToHost);
    cudaMemcpy(S->I_Mag_spin_flip_mp, S_gpu->I_Mag_spin_flip_mp, len_I, cudaMemcpyDeviceToHost);
    cudaMemcpy(S->I_Mag_non_spin_flip_pp, S_gpu->I_Mag_non_spin_flip_pp, len_I, cudaMemcpyDeviceToHost);
    cudaMemcpy(S->I_Mag_non_spin_flip_mm, S_gpu->I_Mag_non_spin_flip_mm, len_I, cudaMemcpyDeviceToHost);
    cudaMemcpy(S->I_Mag_sanspol_p, S_gpu->I_Mag_sanspol_p, len_I, cudaMemcpyDeviceToHost);
    cudaMemcpy(S->I_Mag_sanspol_m, S_gpu->I_Mag_sanspol_m, len_I, cudaMemcpyDeviceToHost);

    // Copy spectral arrays
    cudaMemcpy(S->A_Nuc_unpolarized, S_gpu->A_Nuc_unpolarized, len_A, cudaMemcpyDeviceToHost);
    cudaMemcpy(S->A_Mag_unpolarized, S_gpu->A_Mag_unpolarized, len_A, cudaMemcpyDeviceToHost);
    cudaMemcpy(S->A_Mag_polarized, S_gpu->A_Mag_polarized, len_A, cudaMemcpyDeviceToHost);
    cudaMemcpy(S->A_NucMag, S_gpu->A_NucMag, len_A, cudaMemcpyDeviceToHost);
    cudaMemcpy(S->A_Mag_spin_flip, S_gpu->A_Mag_spin_flip, len_A, cudaMemcpyDeviceToHost);
    cudaMemcpy(S->A_Mag_chiral, S_gpu->A_Mag_chiral, len_A, cudaMemcpyDeviceToHost);
    cudaMemcpy(S->A_Mag_spin_flip_pm, S_gpu->A_Mag_spin_flip_pm, len_A, cudaMemcpyDeviceToHost);
    cudaMemcpy(S->A_Mag_spin_flip_mp, S_gpu->A_Mag_spin_flip_mp, len_A, cudaMemcpyDeviceToHost);
    cudaMemcpy(S->A_Mag_non_spin_flip_pp, S_gpu->A_Mag_non_spin_flip_pp, len_A, cudaMemcpyDeviceToHost);
    cudaMemcpy(S->A_Mag_non_spin_flip_mm, S_gpu->A_Mag_non_spin_flip_mm, len_A, cudaMemcpyDeviceToHost);
    cudaMemcpy(S->A_Mag_sanspol_p, S_gpu->A_Mag_sanspol_p, len_A, cudaMemcpyDeviceToHost);
    cudaMemcpy(S->A_Mag_sanspol_m, S_gpu->A_Mag_sanspol_m, len_A, cudaMemcpyDeviceToHost);

    // === 4. Synchronize to ensure data is ready ===
    cudaDeviceSynchronize();
}

void free_SpectralData(SpectralData* S, SpectralData* S_gpu) {

    free(S->Nq);
    free(S->Ntheta);
    free(S->k_max);
    free(S->dtheta);
    free(S->q);

    free(S->I_Nuc_unpolarized);
    free(S->I_Mag_unpolarized);
    free(S->I_Mag_polarized);
    free(S->I_NucMag);
    free(S->I_Mag_spin_flip);
    free(S->I_Mag_chiral);
    free(S->I_Mag_spin_flip_pm);
    free(S->I_Mag_spin_flip_mp);
    free(S->I_Mag_non_spin_flip_pp);
    free(S->I_Mag_non_spin_flip_mm);
    free(S->I_Mag_sanspol_p);
    free(S->I_Mag_sanspol_m);

    free(S->A_Nuc_unpolarized);
    free(S->A_Mag_unpolarized);
    free(S->A_Mag_polarized);
    free(S->A_NucMag);
    free(S->A_Mag_spin_flip);
    free(S->A_Mag_chiral);
    free(S->A_Mag_spin_flip_pm);
    free(S->A_Mag_spin_flip_mp);
    free(S->A_Mag_non_spin_flip_pp);
    free(S->A_Mag_non_spin_flip_mm);
    free(S->A_Mag_sanspol_p);
    free(S->A_Mag_sanspol_m);

    cudaDeviceSynchronize();
    cudaFree(S_gpu->Nq);
    cudaDeviceSynchronize();
    cudaFree(S_gpu->Ntheta);
    cudaDeviceSynchronize();
    cudaFree(S_gpu->k_max);
    cudaDeviceSynchronize();
    cudaFree(S_gpu->dtheta);
    cudaDeviceSynchronize();
    cudaFree(S_gpu->q);
    cudaDeviceSynchronize();

    cudaFree(S_gpu->I_Nuc_unpolarized);
    cudaDeviceSynchronize();
    cudaFree(S_gpu->I_Mag_unpolarized);
    cudaDeviceSynchronize();
    cudaFree(S_gpu->I_Mag_polarized);
    cudaDeviceSynchronize();
    cudaFree(S_gpu->I_NucMag);
    cudaDeviceSynchronize();
    cudaFree(S_gpu->I_Mag_spin_flip);
    cudaDeviceSynchronize();
    cudaFree(S_gpu->I_Mag_chiral);
    cudaDeviceSynchronize();
    cudaFree(S_gpu->I_Mag_spin_flip_pm);
    cudaDeviceSynchronize();
    cudaFree(S_gpu->I_Mag_spin_flip_mp);
    cudaDeviceSynchronize();
    cudaFree(S_gpu->I_Mag_non_spin_flip_pp);
    cudaDeviceSynchronize();
    cudaFree(S_gpu->I_Mag_non_spin_flip_mm);
    cudaDeviceSynchronize();
    cudaFree(S_gpu->I_Mag_sanspol_p);
    cudaDeviceSynchronize();
    cudaFree(S_gpu->I_Mag_sanspol_m);
    cudaDeviceSynchronize();

    cudaFree(S_gpu->A_Nuc_unpolarized);
    cudaDeviceSynchronize();
    cudaFree(S_gpu->A_Mag_unpolarized);
    cudaDeviceSynchronize();
    cudaFree(S_gpu->A_Mag_polarized);
    cudaDeviceSynchronize();
    cudaFree(S_gpu->A_NucMag);
    cudaDeviceSynchronize();
    cudaFree(S_gpu->A_Mag_spin_flip);
    cudaDeviceSynchronize();
    cudaFree(S_gpu->A_Mag_chiral);
    cudaDeviceSynchronize();
    cudaFree(S_gpu->A_Mag_spin_flip_pm);
    cudaDeviceSynchronize();
    cudaFree(S_gpu->A_Mag_spin_flip_mp);
    cudaDeviceSynchronize();
    cudaFree(S_gpu->A_Mag_non_spin_flip_pp);
    cudaDeviceSynchronize();
    cudaFree(S_gpu->A_Mag_non_spin_flip_mm);
    cudaDeviceSynchronize();
    cudaFree(S_gpu->A_Mag_sanspol_p);
    cudaDeviceSynchronize();
    cudaFree(S_gpu->A_Mag_sanspol_m);
    cudaDeviceSynchronize();
}