#!/bin/bash

# Limpiar archivos de entrada previos
rm -f ../input/p* 

echo "Inicio de la ejecución: $(date)"

for (( idxProton=0; idxProton<=6; idxProton++ ))
do
    for (( idxNeutron=$idxProton; idxNeutron<=6; idxNeutron++ ))
    do
        protonCount=$((idxProton * 2))
        neutronCount=$((idxNeutron * 2))
        echo "Calculando para Np=$protonCount, Nn=$neutronCount"
        
        inputFile="../input/p${protonCount}n${neutronCount}_cin_bigstick.txt"
        rm -f "$inputFile"
        touch "$inputFile"
        
        # Configuración de la entrada para bigstick.x
        echo "e" >> "$inputFile"
        echo "../output/p${protonCount}n${neutronCount}" >> "$inputFile"
        echo "usdbpn" >> "$inputFile"
        echo "$protonCount  $neutronCount" >> "$inputFile"
        echo "0" >> "$inputFile"
        echo "upn" >> "$inputFile"
        echo "usdbpn" >> "$inputFile"
        echo "1 18 $(($protonCount + $neutronCount + 16)) 0.3" >> "$inputFile"
        echo "1 1 1 1 1" >> "$inputFile"
        echo "end" >> "$inputFile"
        echo "ld" >> "$inputFile"
        echo "5 2000" >> "$inputFile"
        
        # Ejecutar bigstick.x y medir el tiempo
        time ./bigstick.x < "$inputFile"
    done
done

Este script produce archivos de salida en la ruta ../output/pXnY.res, donde X e Y representan el número de protones y neutrones, respectivamente. El formato típico del archivo de resultados es el siguiente:

  BIGSTICK Version 7.8.4 Apr  2018
  single-particle file = usdbpn
           4           4
           0 +
         120  iterations

  State      E        Ex         J       T
    1    -87.10445   0.00000    -0.000   0.000
    2    -85.60215   1.50230     2.000   0.000
    3    -82.98830   4.11615     2.000   0.000
    4    -82.73201   4.37243     4.000   0.000
    5    -82.03408   5.07037     3.000   0.000

Para extraer las energías del estado fundamental de todos los núcleos par-par, se implementa una función en C++ que lee los archivos de resultados. Debido a la simetría de isospín en la interacción del modelo de capas, solo se calculan los núcleos con número de neutrones mayor o igual al de protones, y luego se asignan las energías a sus núcleos espejo.

#include <iostream>
#include <fstream>
#include <sstream>
#include <string>

void ObtenerEnergiasFundamentales(double** energiasFondos) {
    const std::string directorio = "../p";
    const std::string extension = ".res";
    double energia;

    for (int protones = 0; protones <= 12; protones += 2) {
        for (int neutrones = protones; neutrones <= 12; neutrones += 2) {
            std::string nombreArchivo = directorio;
            // Construir el nombre del archivo con formato dinámico
            if (protones < 10) nombreArchivo += static_cast<char>(48 + protones);
            else nombreArchivo += std::to_string(protones);
            nombreArchivo += "n";
            if (neutrones < 10) nombreArchivo += static_cast<char>(48 + neutrones);
            else nombreArchivo += std::to_string(neutrones);
            nombreArchivo += extension;

            std::cout << "Procesando: Np=" << protones << ", Nn=" << neutrones << ", Archivo=" << nombreArchivo << std::endl;

            std::ifstream archivo(nombreArchivo);
            std::string linea, palabraClave;
            bool encontrado = false;

            // Buscar la línea que comienza con "State"
            while (std::getline(archivo, linea)) {
                std::stringstream flujo(linea);
                flujo >> palabraClave;
                if (palabraClave == "State") {
                    encontrado = true;
                    break;
                }
            }

            if (encontrado) {
                // Leer la siguiente línea para obtener la energía del primer estado
                std::getline(archivo, linea);
                std::stringstream flujo(linea);
                flujo >> palabraClave >> energia;
                std::cout << "Energía del estado fundamental: " << energia << std::endl;

                // Almacenar en la matriz usando índices basados en la mitad de los números
                int idxP = protones / 2;
                int idxN = neutrones / 2;
                energiasFondos[idxP][idxN] = energia;
                energiasFondos[idxN][idxP] = energia;  // Asignar al núcleo espejo
            }
            archivo.close();
        }
    }
}

Esta función almacena las energías del estado fundamental en una matriz cuadrada energiasFondos, indexada por Np/2 y Nn/2, cubriendo todos los casos par-par dentro de la capa sd. La salida se puede utilizar para análisis posteriores de propiedades nucleares.