fork download
copy 
  1. /*--------------------------------------------------------------
  2.     SIMULACIÓN GALÁCTICA: programa que simula la evolución
  3.     de una galaxia desde unas condiciones iniciales aleatorias
  4.     haciendo uso del algoritmo de velocidades de Verlet.
  5.     La galaxia está formada por muchos sistemas solares
  6.     y por un agujero negro supermasivo central. 
  7.     El objetivo será ver como giran alrededor de dicho agujero
  8.     negro, alcanzando un estado estacionario.
  9. ----------------------------------------------------------------*/
  10.  
  11. #include <iostream>
  12. #include <random>
  13. #include <cmath>
  14. #include <fstream>
  15.  
  16. #define numsist 10 //Número de sistemas solares
  17. #define G 6.351695379e-10 //Constante de grav en kly^3/(M_sol*My^2)
  18. //Donde G=(6.67e-11 m^3/(kgs^2))*(1.989e30 kg/M_sol)*
  19. //*(365.25*24*60*60e6 s/My)^2/(c*365.25*24*60*60e3 m/kly)^3
  20. #define cte_vel 148.6271359 // 1 kly/My= 148627.1359 m/s = 148.6271359 km/s
  21.                             //porque prefiero escribir en mi fichero en km/s
  22.  
  23.  
  24. /*CLASES*/
  25.  
  26.  
  27. //-----------------------------
  28.  
  29.  
  30. //CLASE de los sistemas solares
  31. class cSistema
  32. {
  33. private:
  34.     //Características de los sistemas (comienza por s de Sistema Solar)
  35.     double sPosicionX, sVelocidadX, sPosicionY;
  36.     double  sVelocidadY, sMasa, sExcentricidad;
  37. public:
  38.     cSistema(void); //Constructor
  39.     bool enAgujero(void); //Método para comprobar si ha sido absorbido
  40.     //Métodos para la muestra de los datos de posiciones
  41.     double getPosX(void){return sPosicionX;}; double getPosY(void){return sPosicionY;};
  42.     //Métodos para la muestra de los datos de velocidades
  43.     double getVelX(void){return sVelocidadX;}; double getVelY(void){return sVelocidadY;};
  44.     //Método para la muestra de la masa
  45.     double masa(){return sMasa;};
  46.     //Métodos para establecer las posiciones y las velocidades
  47.     void setPosX(double px){sPosicionX=px;}; void setPosY(double py){sPosicionY=py;};
  48.     void setVelX(double vx){sVelocidadX=vx;}; void setVelY(double vy){sVelocidadY=vy;};
  49.     ~cSistema(void){}; //Destructor
  50. };
  51.  
  52.  
  53. //Función global que me genera un número aleatorio con una semilla diferente
  54. //cada vez que el programa se inicia
  55. std::mt19937& prng() {
  56.     // seed only one PRNG per thread that needs is:
  57.     static thread_local std::mt19937 gen(std::random_device{}());
  58.     return gen;
  59. }
  60.  
  61.  
  62. //MÉTODOS DE cSistema
  63.  
  64. //Constructor de los sistemas solares
  65. cSistema::cSistema(void)
  66. {
  67.     //La galaxia mide unos 200 kly, distibuyo los
  68.     //sistemas en un circulo de radio 100 kly aleatoriamente
  69.     double a=100.0; double pi=3.14159265;
  70.  
  71.     //Distribución de radios desde a/8 hasta a
  72.     std::uniform_real_distribution<double> distribucion(a/8,a);
  73.     double r = distribucion(prng());
  74.     //Distribución de ángulos de 0 a 2pi (circulo completo)
  75.     std::uniform_real_distribution<double> distribucion_phi(0.0,2.0*pi);
  76.     double phi=distribucion_phi(prng());
  77.  
  78.  
  79.     /*Inicialización de parámetros de los sistemas*/
  80.     //La masa de cada sistema solar en masas solares
  81.     sMasa=1.0014;
  82.     //Posiciones aleatorias entre -100 y 100 kly
  83.     sPosicionX= r*cos(phi);
  84.     sPosicionY= r*sin(phi);
  85.     //Velocidades aleatorias entre 0 y cte_vel=0.667 kly/My (=200km/s)
  86.     double rotacion=-1; // 1: clockwise   -1: counter-clockwise 
  87.     double v=sqrt(G*4.154e6/r); //Le doy justo la velocidad orbital, para que trace orbitas círculares
  88.     sVelocidadX=rotacion*v*sin(phi);
  89.     sVelocidadY=-rotacion*v*cos(phi);
  90.  
  91.     //Control posiciones y velocidades por consola
  92.     //std::cout << "Posicion ("<<sPosicionX<<", "<<sPosicionY<<"), velocidad ("<<sVelocidadX<<", "
  93.     //<<sVelocidadY<<")"<<std::endl;
  94.     //La excentricidad se obtiene a partir de las velocidades
  95.     //ver "Vector de excentricidad"
  96.     }
  97.  
  98. //Método para saber si el sistema solar ha sido
  99. //absorbido por el agujero negro (SIN USO, aún)
  100. bool cSistema::enAgujero(void)
  101. {
  102.     if (sPosicionX<2e-5 && sPosicionY<2e-5) //R_aguj=17R_sol=2e-5 klyr
  103.         return true; //Devuelve verdad si está dentro
  104.     else return false; //Falso si no
  105. }
  106.  
  107.  
  108.  
  109. //-----------------------------
  110.  
  111.  
  112. //CLASE del agujero negro
  113. class cAgujero
  114. {
  115. private:
  116.     //Características del agujero negro (comienza por a de Agujero Negro)
  117.     double aPosicionX, aVelocidadX, aPosicionY;
  118.     double  aVelocidadY, aMasa, aExcentricidad;
  119. public:
  120.     cAgujero(void); //Constructor
  121.     //Métodos para la muestra de los datos de posiciones (BH=blackhole)
  122.     double getPosXBH(void){return aPosicionX;}; double getPosYBH(void){return aPosicionY;};
  123.     //Métodos para la muestra de los datos de velocidades (BH=blackhole)
  124.     double getVelXBH(void){return aVelocidadX;}; double getVelYBH(void){return aVelocidadY;};
  125.     //Método para la muestra de la masa (BH=blackhole)
  126.     double masaBH(){return aMasa;};
  127.     //Métodos para establecer las posiciones, las velocidades y la masa
  128.     void setPosXBH(double pxBH){aPosicionX=pxBH;}; void setPosYBH(double pyBH){aPosicionY=pyBH;};
  129.     void setVelXBH(double vxBH){aVelocidadX=vxBH;}; void setVelYBH(double vyBH){aVelocidadY=vyBH;};
  130.     void setMasaBH(double mBH){aMasa=mBH;};
  131.     ~cAgujero(void){}; //Destructor
  132. };
  133.  
  134. //MÉTODOS DE cAgujero
  135.  
  136. //Constructor del agujero negro
  137. cAgujero::cAgujero(void)
  138. {
  139.     /*Valores iniciales para el agujero negro*/
  140.     aMasa=4.154e6; //La masa del agujero negro en masas solares
  141.     aPosicionX=0; aPosicionY=0; //Lo coloco en el centro
  142.     aVelocidadX=0; aVelocidadY=0; //Y sin velocidad
  143. }
  144.  
  145.  
  146. //-----------------------------
  147.  
  148.  
  149.  
  150. //CLASE de la galaxia
  151. class cGalaxia {
  152. private:
  153.     cSistema Solar[numsist];
  154.     cAgujero Negro;
  155.  
  156. public :
  157.    cGalaxia(){};
  158.    ~cGalaxia(){}; //Constructor y destructor
  159.    //Método del algoritmo de verlet
  160.    void verlet(cSistema (&Solar)[numsist], cAgujero Negro,
  161.    double (&ax)[numsist+1], double (&ay)[numsist+1], double h);
  162.    //Método que inicializa las aceleraciones (verlet, pero solo un uso)
  163.    void AceleracionesIniciales(cSistema (&Solar)[numsist], cAgujero Negro,
  164.    double (&ax)[numsist+1], double (&ay)[numsist+1]);
  165.    //HACER FUNCION choques(); //Método de los choques inelésticos
  166. };
  167.  
  168. //Método de las aceleraciones iniciales
  169. //ACELERACION AGUJERO NEGRO ("elemento 1000")
  170. void cGalaxia::AceleracionesIniciales(cSistema (&Solar)[numsist], cAgujero Negro,
  171. double (&ax)[numsist+1], double (&ay)[numsist+1])
  172.  {
  173.     int i, j;
  174.     double xaux,yaux,denom;
  175.     double xauxBH, yauxBH, denomBH;
  176.     //Aceleración de cada sistema
  177.     ax[numsist]=0; ay[numsist]=0; //Inicialización para el agujero negro
  178.     for(i=0;i<numsist;i++) //Evalúo cada sistema
  179.     {
  180.         ax[i]=0;ay[i]=0; //Inicialización aceleraciones 
  181.         for(j=0;j<numsist;j++) //Evaluo la interacción con cada otro sistema
  182.         if (j!=i && (abs(Solar[i].getPosX()-Solar[j].getPosX())<5) 
  183.         && (abs(Solar[i].getPosY()-Solar[j].getPosY())<5)) //Cuadrado de 5x5 para evitar
  184.                                                            //calculos inútiles, solo influencia cercana
  185.         {
  186.             //Variables auxiliares
  187.             xaux=Solar[i].getPosX()-Solar[j].getPosX();
  188.             yaux=Solar[i].getPosY()-Solar[j].getPosY();
  189.             denom=pow(((xaux*xaux)+(yaux*yaux)),1.5);
  190.  
  191.             //Aceleraciones
  192.             ax[i]-=G*Solar[j].masa()*xaux/denom;
  193.             ay[i]-=G*Solar[j].masa()*yaux/denom;
  194.         }
  195.         //Añado la contribución del agujero negro a cada sistema i
  196.         xaux=Solar[i].getPosX()-Negro.getPosXBH();
  197.         yaux=Solar[i].getPosY()-Negro.getPosYBH();
  198.         denom=pow(((xaux*xaux)+(yaux*yaux)),1.5);
  199.         ax[i]-=G*Negro.masaBH()*xaux/denom;
  200.         ay[i]-=G*Negro.masaBH()*yaux/denom;
  201.  
  202.         //control aceleraciones por consola
  203.         //std::cout<<"Aceleracion x"<<i<<": "<<ax[i]<<std::endl;
  204.         //std::cout<<"Aceleracion y"<<i<<": "<<ay[i]<<std::endl;
  205.  
  206.         //Calculo la aceleración del agujero negro
  207.         xauxBH=Negro.getPosXBH()-Solar[i].getPosX();
  208.         yauxBH=Negro.getPosYBH()-Solar[i].getPosY();
  209.         denomBH=pow(((xauxBH*xauxBH)+(yauxBH*yauxBH)),1.5);
  210.         ax[numsist]-=G*Solar[i].masa()*xauxBH/denomBH;
  211.         ay[numsist]-=G*Solar[i].masa()*yauxBH/denomBH;
  212.     }
  213.     //control aceleraciones por consola
  214.     //std::cout<<"Aceleracion xBH: "<<ax[numsist]<<std::endl;
  215.     //std::cout<<"Aceleracion yBH: "<<ay[numsist]<<std::endl;
  216.     return; 
  217.  }
  218.  
  219.  
  220. //Método de aplicación del algoritmo de verlet a la galaxia completa
  221. void cGalaxia::verlet(cSistema (&Solar)[numsist], cAgujero Negro,
  222.  double (&ax)[numsist+1], double (&ay)[numsist+1], double h)
  223. {
  224.     //Declaración de variables
  225.     double wx[numsist+1],wy[numsist+1],xaux,yaux,denom, xauxBH, yauxBH, denomBH;
  226.     int i,j;
  227.  
  228.     //Posiciones, y velocidades auxiliares
  229.     for(i=0; i<numsist; i++)
  230.     {
  231.         //Cambio de r(t) a r(t+h)
  232.         Solar[i].setPosX(Solar[i].getPosX()+(h*Solar[i].getVelX())+((h*h/2.)*ax[i]));
  233.         Solar[i].setPosY(Solar[i].getPosY()+(h*Solar[i].getVelY())+((h*h/2.)*ay[i]));
  234.         //Velocidades auxiliares
  235.         wx[i]=Solar[i].getVelX()+((h/2.)*ax[i]);
  236.         wy[i]=Solar[i].getVelY()+((h/2.)*ay[i]);
  237.     }
  238.     //Para el agujero negro
  239.     Negro.setPosXBH(Negro.getPosXBH()+(h*Negro.getVelXBH())+((h*h/2.)*ax[numsist]));
  240.     Negro.setPosYBH(Negro.getPosYBH()+(h*Negro.getVelYBH())+((h*h/2.)*ay[numsist]));
  241.     wx[numsist]=Negro.getVelXBH()+((h/2.)*ax[numsist]);
  242.     wy[numsist]=Negro.getVelYBH()+((h/2.)*ay[numsist]);
  243.  
  244.  
  245.     //Aceleraciones
  246.     for(i=0; i<numsist;i++) //Evalúo cada sistema
  247.     {
  248.         ax[i]=0;ay[i]=0; //Inicialización aceleraciones
  249.         for(j=0;j<numsist;j++) //Evalúo la interacción con los demás sistemas
  250.         if (j!=i && (abs(Solar[i].getPosX()-Solar[j].getPosX())<5) 
  251.         && (abs(Solar[i].getPosY()-Solar[j].getPosY())<5)) //Cuadrado de 5x5 para no hacer 
  252.                                                            //calculos inútiles, solo influencia cercana
  253.         {
  254.             //Variables auxiliares
  255.             xaux=Solar[i].getPosX()-Solar[j].getPosX();
  256.             yaux=Solar[i].getPosY()-Solar[j].getPosY();
  257.             denom=pow(((xaux*xaux)+(yaux*yaux)),1.5);
  258.  
  259.             //Aceleraciones
  260.             ax[i]-=G*Solar[j].masa()*xaux/denom;
  261.             ay[i]-=G*Solar[j].masa()*yaux/denom;
  262.         }   
  263.         //Añado la contribución del agujero negro a cada sistema i
  264.         xaux=Solar[i].getPosX()-Negro.getPosXBH();
  265.         yaux=Solar[i].getPosY()-Negro.getPosYBH();
  266.         denom=pow(((xaux*xaux)+(yaux*yaux)),1.5);
  267.         ax[i]-=G*Negro.masaBH()*xaux/denom;
  268.         ay[i]-=G*Negro.masaBH()*yaux/denom;
  269.         //Calculo la aceleración del agujero negro
  270.         xauxBH=Negro.getPosXBH()-Solar[i].getPosX();
  271.         yauxBH=Negro.getPosYBH()-Solar[i].getPosY();
  272.         denomBH=pow(((xauxBH*xauxBH)+(yauxBH*yauxBH)),1.5);
  273.         ax[numsist]-=G*Solar[i].masa()*xauxBH/denomBH;
  274.         ay[numsist]-=G*Solar[i].masa()*yauxBH/denomBH;
  275.     }
  276.  
  277.     //Velocidades
  278.     for(i=0;i<numsist;i++)
  279.     {
  280.         Solar[i].setVelX(wx[i]+((h/2.)*ax[i]));
  281.         Solar[i].setVelY(wy[i]+((h/2.)*ay[i]));
  282.     }
  283.     //Para el BH
  284.     Negro.setVelXBH(wx[numsist]+((h/2.)*ax[numsist]));
  285.     Negro.setVelYBH(wy[numsist]+((h/2.)*ay[numsist]));
  286.  
  287.     return;
  288. }
  289.  
  290.  
  291. //=====================================================
  292.  
  293.  
  294.  
  295.  
  296.  
  297. //=====================================================
  298.  
  299.  
  300.  
  301.  
  302.  
  303. /*FUNCIÓN PRINCIPAL*/
  304.  
  305. int main (void)
  306. {
  307.     //DECLARACIÓN DE VARIABLES PRINCIPALES
  308.     cSistema sistemas[numsist]; //Todos los sistemas solares
  309.     cAgujero agujero;
  310.     cGalaxia galaxia;
  311.     int i, j, iter, escribo;
  312.     double aceleracionx[numsist+1], aceleraciony[numsist+1];
  313.     double t_fin,t_aux, h;
  314.     std::ofstream pos; //Llamo a mi fichero de posiciones "pos"
  315.     std::ofstream vel; //Llamo a mi fichero de velocidades "vel"
  316.  
  317.     //Puesta a punto de los ficheros
  318.     pos.open("posiciones.txt"); //Lo abro en un fichero llamado posiciones.txt
  319.     if (!pos.is_open()) std::cout<<"Error al abrir el fichero posiciones.txt";
  320.     vel.open("velocidades.txt"); //Lo abro en un fichero llamado velocidades.txt
  321.     if (!vel.is_open()) std::cout<<"Error al abrir el fichero velocidades.txt";
  322.     pos << "#PosicionX   PosicionY"<<std::endl; //Almohadilla porque es comentario en gnuplot
  323.     vel << "#Radio    Velocidad"<<std::endl; //Almohadilla porque es comentario en gnuplot
  324.  
  325.     //Lectura de datos temporales por pantalla
  326.     /*std::cout<<"Introduzca el paso infinitesimal temporal: ";
  327.     std::cin>>h;std::cout<<std::endl;
  328.     std::cout<<"Introduzca el tiempo de evolucion galactica: ";
  329.     std::cin>>t_fin;std::cout<<std::endl;
  330.     std::cout<<"Escribir datos en ficheros cada cuantas iteraciones: ";
  331.     std::cin>>escribo;std::cout<<std::endl;*/
  332.  
  333.     //Por si me cansa la escritura en la consola:
  334.     //Año galactico ~250My
  335.     h=0.005;t_fin=2e4;
  336.     escribo=(int) t_fin/(h*1e3);
  337.  
  338.     //Implementación del algoritmo a la función principal, INICIALIZACIÓN
  339.     galaxia.AceleracionesIniciales(sistemas,agujero,aceleracionx,aceleraciony);
  340.  
  341.     //Para simplificar las cosas, uso numeros enteros
  342.     //Porque me es útil para saber el porcentaje de ejecución
  343.     //y tambien el numero de escrituras en archivo (va con division entera)
  344.     int hmod=lround(h/h);
  345.     int t_auxmod=lround(t_aux/h); 
  346.     int t_finmod=lround(t_fin/h);
  347.     int orden=floor(log10(t_fin));
  348.  
  349.     for(t_auxmod=hmod; t_auxmod<=t_finmod+hmod; t_auxmod=t_auxmod+hmod) //de h a t_fin con pasos de h
  350.     {
  351.         galaxia.verlet(sistemas,agujero,aceleracionx,aceleracionx,h); //ALGORITMO
  352.  
  353.         //ESCRITURA
  354.         //El numero de iteraciones del algoritmo será:
  355.         iter=t_auxmod/hmod;
  356.         //Escribimos cada determinado número de iteraciones
  357.         if((iter%escribo)==0){
  358.         //Escritura en fichero pos y vel del agujero negro
  359.         pos << agujero.getPosXBH() << "\t" << agujero.getPosYBH() << std::endl;
  360.             vel << sqrt(agujero.getPosXBH()*agujero.getPosXBH()+agujero.getPosYBH()*agujero.getPosYBH())
  361.         << "\t" <<sqrt(agujero.getVelXBH()*agujero.getVelXBH()+agujero.getVelYBH()*agujero.getVelYBH())*cte_vel<<std::endl;
  362.         //Escritura de los sistemas solares
  363.         for(i=0;i<numsist;i++){
  364.         pos<<sistemas[i].getPosX() <<"\t"<< sistemas[i].getPosY() << std::endl;  
  365.             vel<<sqrt(sistemas[i].getPosX()*sistemas[i].getPosX()+sistemas[i].getPosY()*sistemas[i].getPosY())
  366.         <<"\t"<<sqrt(sistemas[i].getVelX()*sistemas[i].getVelX()+sistemas[i].getVelY()*sistemas[i].getVelY())*cte_vel<<std::endl;
  367.         }//endfor
  368.         pos<<std::endl<<std::endl;vel<<std::endl<<std::endl;//Dos espacios en blanco separan cada iteración en ambos ficheros
  369.         }//endif
  370.  
  371.  
  372.         //PORCENTAJE DE EJECUCIÓN
  373.         int pt=10; //Paso porcentual que se muestra en pantalla
  374.         //Exponentes
  375.         int grande=1+orden;
  376.         int exp1=2+grande; int exp2=1+grande; int exp3=grande;
  377.         //Porcentaje de ejecución, cada 10%, fallos por ser numeros inexactos y redondeo de int
  378.         unsigned long int porcentaje=100*t_aux/t_fin;
  379.         if(porcentaje%pt==0){ //Detección paso a paso no muy precisa
  380.         porcentaje = pow(pt,exp1)*t_auxmod/t_finmod; //Para ser más preciso en mis cálculos (division de enteros redondea)
  381.         if(porcentaje%(int)pow(pt,exp2)==0) std::cout << "Ejecutado al "<<porcentaje/pow(pt,exp3)<<"%"<<std::endl;
  382.         //^Detección precisa
  383.         }//endif
  384.     }
  385.  
  386.     //Cierro los ficheros
  387.     pos.close();
  388.     vel.close();
  389.  
  390.     return 0;
  391. }
  392.  
Time limit exceeded #stdin #stdout 5s 4232KB
comments (?)
stdin
copy 
Standard input is empty
stdout
copy 
Error al abrir el fichero posiciones.txtError al abrir el fichero velocidades.txtEjecutado al 10%
Ejecutado al 20%
Ejecutado al 30%
Ejecutado al 40%
Ejecutado al 50%
Ejecutado al 60%
Ejecutado al 70%
Ejecutado al 80%
https://ideone.com/8awfhm
language:
C++ (gcc 8.3)
created:
3 years ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!