KG-ICN.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d66bb16d-b338-476d-8ceb-199bbc7c9227",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Importamos bibliotecas.\n",
    "using Plots, LaTeXStrings, DelimitedFiles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0f354e7c-3e54-4c2d-bb84-fbd3986b914a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculamos la condición inicial de Carlos.\n",
    "\n",
    "# Este es el sistema estático de Carlos.\n",
    "function rhsStatic(w, r, a, α, ϕ, ψ)\n",
    "    \n",
    "    # Definimos las ecuaciones diferenciales de los estados estáticos.\n",
    "    da = (1 - a^2)/(2*r) + (r/2)*(ψ^2 + (a*ϕ)^2*(1 + (w/α)^2))\n",
    "    dα = (a^2 - 1)/r + da - r*(a*ϕ)^2\n",
    "    dϕ = ψ\n",
    "    dψ = -(2/r + dα - da)*ψ + a^2*(1 - (w/α)^2)*ϕ\n",
    "\n",
    "    return [a*da, α*dα, dϕ, dψ]\n",
    "end\n",
    "\n",
    "# Definimos una función que haga rk4 para calcular el siguiente paso.\n",
    "function rk4Static(dr, w, r, a, α, ϕ, ψ)\n",
    "\n",
    "    k1  = rhsStatic(w, r, a, α, ϕ, ψ)\n",
    "    k1a = k1[1]\n",
    "    k1α = k1[2]\n",
    "    k1ϕ = k1[3]\n",
    "    k1ψ = k1[4]\n",
    "\n",
    "    k2  = rhsStatic(w, r + dr/2, a + (k1a*dr/2), α + (k1α*dr/2), ϕ + (k1ϕ*dr/2), ψ + (k1ψ*dr/2))\n",
    "    k2a = k2[1]\n",
    "    k2α = k2[2]\n",
    "    k2ϕ = k2[3]\n",
    "    k2ψ = k2[4]\n",
    "\n",
    "    k3  = rhsStatic(w, r + dr/2, a + (k2a*dr/2), α + (k2α*dr/2), ϕ + (k2ϕ*dr/2), ψ + (k2ψ*dr/2))\n",
    "    k3a = k3[1]\n",
    "    k3α = k3[2]\n",
    "    k3ϕ = k3[3]\n",
    "    k3ψ = k3[4]\n",
    "\n",
    "    k4  = rhsStatic(w, r + dr, a + k3a*dr, α + k3α*dr, ϕ + k3ϕ*dr, ψ + k3ψ*dr)\n",
    "    k4a = k4[1]\n",
    "    k4α = k4[2]\n",
    "    k4ϕ = k4[3]\n",
    "    k4ψ = k4[4]\n",
    "\n",
    "    rrk4 = r + dr\n",
    "    ark4 = a + (k1a + 2*k2a + 2*k3a + k4a)*(dr/6)\n",
    "    αrk4 = α + (k1α + 2*k2α + 2*k3α + k4α)*(dr/6)\n",
    "    ϕrk4 = ϕ + (k1ϕ + 2*k2ϕ + 2*k3ϕ + k4ϕ)*(dr/6)\n",
    "    ψrk4 = ψ + (k1ψ + 2*k2ψ + 2*k3ψ + k4ψ)*(dr/6)\n",
    "    \n",
    "    return [rrk4, ark4, αrk4, ϕrk4, ψrk4]\n",
    "end\n",
    "\n",
    "function solver(ϕi, w, Nr, dr)\n",
    "    \n",
    "    # Definimos las condiciones iniciales.\n",
    "    r0 = dr\n",
    "    a0 = 1\n",
    "    α0 = 1\n",
    "    ϕ0 = ϕi\n",
    "    ψ0 = 0.0\n",
    "\n",
    "    # Definimos los arreglos a usar.\n",
    "    r = zeros(Nr)\n",
    "    a = zeros(Nr)\n",
    "    α = zeros(Nr)\n",
    "    ϕ = zeros(Nr)\n",
    "    ψ = zeros(Nr)\n",
    "\n",
    "    # Llenamos las condiciones iniciales (y las simetrías).\n",
    "    r[2] = r0\n",
    "    a[2] = a0\n",
    "    α[2] = α0\n",
    "    ϕ[2] = ϕ0\n",
    "    ψ[2] = ψ0\n",
    "\n",
    "    r[1] = -r0\n",
    "    a[1] = +a0\n",
    "    α[1] = +α0\n",
    "    ϕ[1] = +ϕ0\n",
    "    ψ[1] = -ψ0\n",
    "\n",
    "    # Comenzamos a solucionar en sí.\n",
    "    for j in 2:Nr-1\n",
    "        soluto = rk4Static(dr, w, r[j], a[j], α[j], ϕ[j], ψ[j])\n",
    "\n",
    "        r[j+1] = soluto[1]\n",
    "        a[j+1] = soluto[2]\n",
    "        α[j+1] = soluto[3]\n",
    "        ϕ[j+1] = soluto[4]\n",
    "        ψ[j+1] = soluto[5]\n",
    "    end\n",
    "\n",
    "    return [r, a, α, ϕ, ψ]\n",
    "end\n",
    "\n",
    "function shooting(ϕi, w, Nr, dr)\n",
    "    \n",
    "    # Definimos la primer diferencial de 'w'.\n",
    "    tolerancia = dr^4\n",
    "    dw = 0.000001\n",
    "\n",
    "    # Definimos un contador para ver cuántas veces necesitamos solucionar.\n",
    "    counter = 0\n",
    "    \n",
    "    # Solucionamos por primera vez.\n",
    "    ϕ = solver(ϕi, w, Nr, dr)[4]\n",
    "    counter += 1\n",
    "\n",
    "    while abs( ϕ[end] ) > tolerancia\n",
    "        \n",
    "        if ϕ[end] > tolerancia\n",
    "            # Actualizamos la 'w'.\n",
    "            w = w + dw\n",
    "\n",
    "            # Volvemos a solucionar.\n",
    "            ϕ = solver(ϕi, w, Nr, dr)[4]\n",
    "\n",
    "        else\n",
    "            # Aquí la 'w' ha sido muy grande, entonces debemos regresar.\n",
    "            w = w - dw\n",
    "\n",
    "            # Ahora hacemos más pequeña la varaición y volvemos a avanzar.\n",
    "            dw = dw/2\n",
    "            w = w + dw\n",
    "\n",
    "            # Volvemos a solucionar.\n",
    "            ϕ = solver(ϕi, w, Nr, dr)[4]\n",
    "        end\n",
    "        \n",
    "        counter += 1\n",
    "    end\n",
    "        \n",
    "    # Solucionamos una última vez (pero no la contamos).\n",
    "    solucion = solver(ϕi, w, Nr, dr)\n",
    "    r = solucion[1]\n",
    "    a = solucion[2]\n",
    "    α = solucion[3]\n",
    "    ϕ = solucion[4]\n",
    "    ψ = solucion[5]\n",
    "    \n",
    "    return [counter, w, r, a, α, ϕ, ψ]\n",
    "end\n",
    "\n",
    "function initial(ϕi, w, Nr, dr)\n",
    "\n",
    "    # Aplicamos el shooting.\n",
    "    soluto = shooting(ϕi, w, Nr, dr)\n",
    "    counter = soluto[1]\n",
    "    w = soluto[2]\n",
    "    r = soluto[3]\n",
    "    a = soluto[4]\n",
    "    α = soluto[5]\n",
    "    ϕ = soluto[6]\n",
    "    ψ = soluto[7]\n",
    "    println(\"Iteraciones del shooting: $counter \\n\")\n",
    "    \n",
    "    # Redimensionamos alpha.\n",
    "    α00 = 1/(a[Nr]*α[Nr])\n",
    "    α = α00 .* α\n",
    "\n",
    "    # Calculamos la parte real e imaginaria de las funciones auxiliares.\n",
    "    ϕ1 = copy(ϕ)\n",
    "    ϕ2 = zeros(Nr)\n",
    "\n",
    "    ψ1 = copy(ψ)\n",
    "    ψ2 = zeros(Nr)\n",
    "\n",
    "    π1 = zeros(Nr)\n",
    "    π2 = -(w*α00) .* ϕ .* ( a ./ α )\n",
    "\n",
    "    m = (r ./ 2) .* ( 1 .- (1 ./ a .^ 2) )\n",
    "    σ = a .* α\n",
    "    \n",
    "    return [r, m, a, α, σ, ϕ1, ϕ2, ψ1, ψ2, π1, π2]\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7c8411c4-33af-48f4-94b0-8c588b5bf517",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Definimos una función para los rhs.\n",
    "function rhs(Nr, dr, r, a, α, ϕ1, ϕ2, ψ1, ψ2, π1, π2)\n",
    "\n",
    "    # Definimos arreglos para guardar los lados derechos (los rhs: right-hand side).\n",
    "    ϕ1_f = zeros(Nr)\n",
    "    ϕ2_f = zeros(Nr)\n",
    "    \n",
    "    ψ1_f = zeros(Nr)\n",
    "    ψ2_f = zeros(Nr)\n",
    "\n",
    "    π1_f = zeros(Nr)\n",
    "    π2_f = zeros(Nr)\n",
    "\n",
    "    # Llenamos estos arreglos con (9) de Argelia-Siddhartha.\n",
    "    for j in 1:Nr\n",
    "        ϕ1_f[j] = (α[j]/a[j])*π1[j]\n",
    "        ϕ2_f[j] = (α[j]/a[j])*π2[j]\n",
    "    end\n",
    "\n",
    "    for j in 2:Nr-1\n",
    "        ψ1_f[j] = ( (α[j+1]/a[j+1])*π1[j+1] - (α[j-1]/a[j-1])*π1[j-1] )/(2*dr)\n",
    "        ψ2_f[j] = ( (α[j+1]/a[j+1])*π2[j+1] - (α[j-1]/a[j-1])*π2[j-1] )/(2*dr)\n",
    "    end\n",
    "\n",
    "    for j in 2:Nr-1\n",
    "        π1_f[j] = (1/r[j]^2)*( r[j+1]^2*(α[j+1]/a[j+1])*ψ1[j+1] -\n",
    "                    r[j-1]^2*(α[j-1]/a[j-1])*ψ1[j-1] )/(2*dr) -\n",
    "                    a[j]*α[j]*ϕ1[j]\n",
    "\n",
    "        π2_f[j] = (1/r[j]^2)*( r[j+1]^2*(α[j+1]/a[j+1])*ψ2[j+1] -\n",
    "                    r[j-1]^2*(α[j-1]/a[j-1])*ψ2[j-1] )/(2*dr) -\n",
    "                    a[j]*α[j]*ϕ2[j]\n",
    "    end\n",
    "\n",
    "    return [ϕ1_f, ϕ2_f, ψ1_f, ψ2_f, π1_f, π2_f]\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0d945865-6569-4a83-a918-46071af20d60",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Definimos una función para solucionar la métrica.\n",
    "\n",
    "# Las ecuaciones son (10) y (11) de Argelia-Siddhartha.\n",
    "function metricrk2(Nr, dr, r, ϕ1, ϕ2, ψ1, ψ2, π1, π2)\n",
    "\n",
    "    #------------------------------------------------\n",
    "    # Solucionamos primero para a.\n",
    "    a = zeros(Nr)\n",
    "    a[2] = 1\n",
    "\n",
    "    for j in 2:Nr-1\n",
    "        \n",
    "        rho = (r[j]/2)*( π1[j]^2 + π2[j]^2 + ψ1[j]^2 + ψ2[j]^2 +\n",
    "                (a[j]^2)*(ϕ1[j]^2 + ϕ2[j]^2) )\n",
    "        \n",
    "        geo = ( a[j]^2 - 1 )/( 2*r[j] )\n",
    "\n",
    "        aux_p = rho - geo\n",
    "\n",
    "        k1  = a[j]*aux_p\n",
    "\n",
    "        x   = r[j] + dr/2\n",
    "        y   = a[j] + k1*dr/2\n",
    "\n",
    "        #------------------------------------------------\n",
    "\n",
    "        rho = (x/2)*( π1[j+1]^2 + π2[j+1]^2 + ψ1[j+1]^2 + ψ2[j+1]^2 +\n",
    "                y^2*(ϕ1[j+1]^2 + ϕ2[j+1]^2) )\n",
    "        \n",
    "        geo = ( y^2 - 1 )/( 2*x )\n",
    "\n",
    "        aux = rho - geo\n",
    "\n",
    "        k2  = y*( aux + aux_p )/2\n",
    "\n",
    "        #------------------------------------------------\n",
    "        a[j+1] = a[j] + k2*dr\n",
    "\n",
    "    end\n",
    "\n",
    "    a[1] = a[2]\n",
    "\n",
    "    #------------------------------------------------\n",
    "    # Ahora solucionamos para alpha.\n",
    "    # Necesitamos solucionar hacia atrás porque conocemos una condición\n",
    "    # en la frontera exterior, no en el inicio.\n",
    "    α = zeros(Nr)\n",
    "    α[Nr] = 1/a[Nr]\n",
    "\n",
    "    for j in Nr-1:-1:2\n",
    "\n",
    "        rho = (r[j+1]/2)*( π1[j+1]^2 + π2[j+1]^2 + ψ1[j+1]^2 + ψ2[j+1]^2 -\n",
    "                a[j+1]^2*(ϕ1[j+1]^2 + ϕ2[j+1]^2) )\n",
    "\n",
    "        geo = ( a[j+1]^2 - 1 )/( 2*r[j+1] )\n",
    "\n",
    "        aux_p = geo + rho\n",
    "\n",
    "        k1 = α[j+1]*aux_p\n",
    "\n",
    "        x  = r[j+1] - dr/2\n",
    "        y  = α[j+1] - k1*dr/2\n",
    "\n",
    "        #------------------------------------------------\n",
    "        rho = (x/2)*( π1[j]^2 + π2[j]^2 + ψ1[j]^2 + ψ2[j]^2 -\n",
    "                a[j]^2*(ϕ1[j]^2 + ϕ2[j]^2) )\n",
    "\n",
    "        geo = ( a[j]^2 - 1 )/( 2*x )\n",
    "\n",
    "        aux = geo + rho\n",
    "\n",
    "        k2  = y*( aux + aux_p )/2\n",
    "\n",
    "        #------------------------------------------------\n",
    "        α[j] = α[j+1] - k2*dr\n",
    "        \n",
    "    end\n",
    "\n",
    "    α[1] = α[2]\n",
    "    \n",
    "    return [a, α]\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d63d5261-6276-42bc-b528-45af683db8d1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Definimos una función para el ICN en sí.\n",
    "function icn(Nr, dr, dt, r, a, α, ϕ1, ϕ2, ψ1, ψ2, π1, π2)\n",
    "\n",
    "    # Limpiamos adot.\n",
    "    adot = zeros(Nr)\n",
    "\n",
    "    # Guardamos la información previa.\n",
    "    a_p  = copy(a)\n",
    "    α_p  = copy(α)\n",
    "    ϕ1_p = copy(ϕ1)\n",
    "    ϕ2_p = copy(ϕ2)\n",
    "    ψ1_p = copy(ψ1)\n",
    "    ψ2_p = copy(ψ2)\n",
    "    π1_p = copy(π1)\n",
    "    π2_p = copy(π2)\n",
    "    \n",
    "    #------------------------------------------------\n",
    "    # Iniciamos el ciclo interno del ICN.\n",
    "    for k in 1:3\n",
    "\n",
    "        # Calculamos el tamaño del paso temporal en función de la iteración.\n",
    "        if k < 3\n",
    "            dtw = dt/2\n",
    "        else\n",
    "            dtw = dt\n",
    "        end\n",
    "\n",
    "        # Calculamos los lados derechos.\n",
    "        fuentes = rhs(Nr, dr, r, a, α, ϕ1, ϕ2, ψ1, ψ2, π1, π2)\n",
    "        ϕ1_f = fuentes[1]\n",
    "        ϕ2_f = fuentes[2]\n",
    "        ψ1_f = fuentes[3]\n",
    "        ψ2_f = fuentes[4]\n",
    "        π1_f = fuentes[5]\n",
    "        π2_f = fuentes[6]\n",
    "\n",
    "        # Actualizamos los campos.\n",
    "        for j in 2:Nr-1\n",
    "            ϕ1[j] = ϕ1_p[j] + dtw*ϕ1_f[j]\n",
    "            ϕ2[j] = ϕ2_p[j] + dtw*ϕ2_f[j]\n",
    "            ψ1[j] = ψ1_p[j] + dtw*ψ1_f[j]\n",
    "            ψ2[j] = ψ2_p[j] + dtw*ψ2_f[j]\n",
    "            π1[j] = π1_p[j] + dtw*π1_f[j]\n",
    "            π2[j] = π2_p[j] + dtw*π2_f[j]\n",
    "        end\n",
    "        \n",
    "        # Imponemos las condiciones en la frontera interior.\n",
    "        ϕ1[1] = ϕ1_p[1] + dtw*ϕ1_f[1]\n",
    "        ϕ2[1] = ϕ2_p[1] + dtw*ϕ2_f[1]\n",
    "        ψ1[1] = -ψ1[2]\n",
    "        ψ2[1] = -ψ2[2]\n",
    "        π1[1] = +π1[2]\n",
    "        π2[1] = +π2[2]\n",
    "\n",
    "        # Imponemos las condiciones en la frontera exterior.\n",
    "        div = 1 + dtw/dr + dtw/r[Nr]\n",
    "        ϕ1[Nr] = ϕ1_p[Nr] + dtw*ϕ1_f[Nr]\n",
    "        ϕ2[Nr] = ϕ2_p[Nr] + dtw*ϕ2_f[Nr]\n",
    "        π1[Nr] = (π1[Nr-1]*(dtw/dr - 1) + π1_p[Nr]*(1 - dtw/dr - dtw/r[Nr]) + π1_p[Nr-1]*(dtw/dr + 1))/div\n",
    "        π2[Nr] = (π2[Nr-1]*(dtw/dr - 1) + π2_p[Nr]*(1 - dtw/dr - dtw/r[Nr]) + π2_p[Nr-1]*(dtw/dr + 1))/div\n",
    "        ψ1[Nr] = -π1[Nr] - ϕ1[Nr]/r[Nr]\n",
    "        ψ2[Nr] = -π2[Nr] - ϕ2[Nr]/r[Nr]\n",
    "\n",
    "        # Calculamos la métrica.\n",
    "        metricas = metricrk2(Nr, dr, r, ϕ1, ϕ2, ψ1, ψ2, π1, π2)\n",
    "        a = metricas[1]\n",
    "        α = metricas[2]        \n",
    "\n",
    "    #------------------------------------------------\n",
    "    # Terminamos el ciclo interno.\n",
    "    end\n",
    "    \n",
    "    # Calculamos la masa y la métrica.\n",
    "    m = (r ./ 2) .* ( 1 .- (1 ./ a .^2) )\n",
    "    σ = a .* α\n",
    "\n",
    "    # Calculamos adot.\n",
    "    Srr_p = r .* α_p .* (ψ1_p .* π1_p .+ ψ2_p .* π2_p)\n",
    "    Srr = r .* α .* (ψ1 .* π1 .+ ψ2 .* π2)\n",
    "\n",
    "    adot = ((a .- a_p) ./ dt) .- (1/2) .* (Srr_p .+ Srr)\n",
    "    \n",
    "    return [m, a, α, σ, ϕ1, ϕ2, ψ1, ψ2, π1, π2, adot]\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b46d125a-d9a1-4678-b763-68770ba6d8e6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Definimos la función main.\n",
    "function main(ϕi, w, rf, dr, tf, courant)\n",
    "    \n",
    "    # Definimos la longitud de los arreglos.\n",
    "    Nr = round(Int, rf/dr)\n",
    "    Nt = round(Int, tf*(Nr/(courant*rf)))\n",
    "    dt = courant*dr\n",
    "\n",
    "    # Calculamos la condición de guardado temporal y espacial.\n",
    "    savedataT = round(Int, Nt/100)\n",
    "    savedataR = round(Int, Nr/100)\n",
    "\n",
    "    # Definimos arreglos que guardarán la información.\n",
    "    adot_dat = []\n",
    "    phi1_dat = []\n",
    "\n",
    "    # Iniciamos el tiempo en cero.\n",
    "    t = 0\n",
    "    \n",
    "    # Calculamos las condiciones iniciales.\n",
    "    iniciales = initial(ϕi, w, Nr, dr)\n",
    "    r  = iniciales[1]\n",
    "    m  = iniciales[2]\n",
    "    a  = iniciales[3]\n",
    "    α  = iniciales[4]\n",
    "    σ  = iniciales[5]\n",
    "    ϕ1 = iniciales[6]\n",
    "    ϕ2 = iniciales[7]\n",
    "    ψ1 = iniciales[8]\n",
    "    ψ2 = iniciales[9]\n",
    "    π1 = iniciales[10]\n",
    "    π2 = iniciales[11]\n",
    "    adot = zeros(Nr)\n",
    "\n",
    "    # Guardamos las condiciones iniciales.\n",
    "    for j in 2:Nr\n",
    "        push!(adot_dat, [t, r[j], adot[j]])\n",
    "    end\n",
    "    push!(phi1_dat, [t, ϕ1[2]])\n",
    "\n",
    "    # Graficamos las condiciones iniciales.\n",
    "    grafica = plot(layout = (1, 3), size = (1500, 500))\n",
    "\n",
    "    plot!(grafica[1], r, adot, label = \"adot(r)\", linecolor=:teal)\n",
    "    \n",
    "    plot!(grafica[2], r, a, label = \"a(r)\", linecolor=:magenta)\n",
    "    plot!(grafica[2], r, σ, label = \"σ(r)\", linecolor=:sienna)\n",
    "    plot!(grafica[2], r, α, label = \"α(r)\", linecolor=:teal)\n",
    "\n",
    "    plot!(grafica[3], r, ϕ1, label = \"ϕ1(r)\", linecolor=:magenta)\n",
    "    plot!(grafica[3], r, ϕ2, label = \"ϕ2(r)\", linecolor=:teal)\n",
    "\n",
    "    # Iniciamos el ciclo principal de evolución.\n",
    "    for n in 2:Nt\n",
    "\n",
    "        # Avanzamos en el tiempo.\n",
    "        t = t + dt\n",
    "\n",
    "        # Hacemos el ICN.\n",
    "        crank = icn(Nr, dr, dt, r, a, α, ϕ1, ϕ2, ψ1, ψ2, π1, π2)\n",
    "        m  = crank[1]\n",
    "        a  = crank[2]\n",
    "        α  = crank[3]\n",
    "        σ  = crank[4]\n",
    "        ϕ1 = crank[5]\n",
    "        ϕ2 = crank[6]\n",
    "        ψ1 = crank[7]\n",
    "        ψ2 = crank[8]\n",
    "        π1 = crank[9]\n",
    "        π2 = crank[10]\n",
    "        adot = crank[11]\n",
    "\n",
    "        # Graficamos y gurdamos.\n",
    "        if n % savedataT == 0\n",
    "            for j in 2:Nr\n",
    "                push!(adot_dat, [t, r[j], adot[j]])\n",
    "            end\n",
    "            push!(phi1_dat, [t, ϕ1[2]])\n",
    "            \n",
    "            plot!(grafica[1], r, adot, label=\"\", linecolor=:teal)\n",
    "            \n",
    "            plot!(grafica[2], r, a, label=\"\", linecolor=:magenta)\n",
    "            plot!(grafica[2], r, σ, label=\"\", linecolor=:sienna)\n",
    "            plot!(grafica[2], r, α, label=\"\", linecolor=:teal)\n",
    "\n",
    "            plot!(grafica[3], r, ϕ1, label=\"\", linecolor=:magenta)\n",
    "            plot!(grafica[3], r, ϕ2, label=\"\", linecolor=:teal)\n",
    "        end\n",
    "\n",
    "    end\n",
    "\n",
    "    # Mostramos la gráfica.\n",
    "    display(grafica)\n",
    "\n",
    "    # Escribimos los archivos en sí.\n",
    "    writedlm(\"adot.dat\", adot_dat)\n",
    "    writedlm(\"phi1.dat\", phi1_dat)\n",
    "    \n",
    "    return\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eddf086e-e6de-4ce1-aa41-c1d0399540e5",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Ejecutamos la función principal.\n",
    "@time begin\n",
    "    ϕi = 0.110             # Parámetro inicial del campo.\n",
    "    w  = 1.0809049555648   # Frecuencia de oscilación.\n",
    "    rf = 45                # Radio de integración final.\n",
    "    dr = 0.01              # Tamaño del paso.\n",
    "    tf = 10                # Tiempo final de integración.\n",
    "    courant = 0.5          # parámetro de courant.\n",
    "        \n",
    "    main(ϕi, w, rf, dr, tf, courant)\n",
    "end"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "67fadf7f-84d8-4f2f-9be1-0817221b7df5",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.10.3",
   "language": "julia",
   "name": "julia-1.10"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.10.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}