/******************************************************************/
/*                                                                */
/* 1-D ACCRETION DISC EVOLUTION MODEL                             */
/*                                                                */
/* Richard Alexander - March 2009; this version May 2011          */
/* [History: 2008 - cleaned up and documented.                    */
/* Original code dates from September 2004; heavily modified      */
/* in 2007/08.]                                                   */  
/*                                                                */
/* INTEGRATES FOLLOWING EQUATION FORWARD IN TIME:                 */
/*                                                                */
/* (1)  dSigma/dt = (3/R)d/dR[sqrt(R)d/dR(nu Sigma sqrt(R))]      */
/*                                                                */
/*                                                                */
/* CODE UNITS ARE Msun, AU & yr   (G = 4pi^2)                     */
/* THIS IS ONLY NEEDED IF SOLVING THE ENERGY EQUATION THOUGH -    *
/* THIS VERSION, WITH CONSTANT nu, IS TECHNICALLY DIMENSIONLESS   */
/*                                                                */
/* GRID IS EQUISPACED IN R^(1/2)                                  */
/* Eq (1) SOLVED USING FIRST-ORDER EXPLCIT METHOD                 */
/*                                                                */
/* VARIABLE CHANGES:                                              */
/*        X = 2*sqrt(R)                                           */
/*        Y = 3*nu*Sigma*sqrt(R)                                  */
/*        S = Sigma*R^(3/2)                                       */
/*                                                                */
/*  => Eq (1) becomes   dS/dt = d2Y/dX2                           */
/*                                                                */
/* AT EACH TIMESTEP, QUANTITIES F[*][0] ARE INTEGRATED FORWARD IN */
/* TIME TO F[*][1].  [1] VALUES ARE THEN OVER-WRITTEN TO THE [0]  */
/* ARRAY AT THE END OF THE LOOP.                                  */ 
/* TIME-STEPPING IS ADAPTIVE, USING A COURANT-LIKE CONDITION      */
/*                                                                */
/* OUTPUT IS A SERIES OF FILES sigma_XXXX.dat, WITH FORMAT:       */
/*    time                                                        */
/*    R[in]      sigma[in]                                        */
/*    R[in+1]    sigma[in+1]                                      */
/*    R[in+2]    sigma[in+2]                                      */
/*    R[in+3]    sigma[in+3]                                      */
/*    ...       ...                                               */
/*    R[out]     sigma[out]                                       */
/*                                                                */
/* FILES CAN EASILY BE READ IN TO ANY STANDARD PLOTTING PACKAGE   */
/*                                                                */
/* ALSO GENERATES A SINGLE FILE mdot_???.dat, WITH GLOBAL INFO AT */
/* EACH OUTPUT STEP.  EACH LINE IS FORMATTED:                     */
/*   time   Mdot    M_disc                                        */
/*                                                                */
/*                                                                */
/* THIS VERSION USES A CONSTANT VISCOSITY AND A "SPIKE" AT R=1 AS */
/* THE INITIAL CONDITIONS - IT'S EFFECTIVELY THE SPREADING RING   */
/* TEST FROM PRINGLE (1981).  THE VISCOSITY LAW AND INITIAL DISC  */
/* PROFILE CAN BE MODIFIED TO ANY ARBITRARY FORM AT LINE 150.     */
/*                                                                */
/******************************************************************/

/* STANDARD LIBRARY FILES */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>

/* MAIN PROGRAM */
int main()
{
  int i,j,k,t,in,out,dump_counter;
  /* COMPUTATIONAL PARAMS */
  static const double epsilon=1E-10;
  static const double Courant=0.05;
  /* GRID & VISC */
  double Rz[10001],X[10001],nu[10001],Rf[10001],XX[10001];
  /* GAS */
  double Y[10001][2],sigma[10001][2],S[10001][2];
  /* OTHERS */
  double M_dot[2],M_acc,M_out;
  double dX,dt,minstep,gstep;
  double disc_mass;
  double time,time_to_dump;
  double dR,nu_s;
  /* STRINGS AND FILE POINTERS */
  char sigfile[100],mfile[100],root[100];
  FILE *sigmaout,*mdot;

  
/******************************************************************/
/*                                                                */
  
  /* USER-DEFINED PARAMETERS              */
  
  /* OUTER BOUNDARY CONDITION:            */
  /* 0 = ZERO-TORQUE ; 1 = REFLECTING     */
  static const int BC_flag=0;
  /* DATA DUMP FREQUENCY  */
  static const double step=1.0E-3;
  /* STOPPING TIME FOR INTEGRATION  */
  static const double t_max=10.0;
  
  /* BASE STRING FOR MAKING OUTPUT FILENAMES */
  sprintf(root,"1dring");

/*                                                                */
/*                                                                */
/******************************************************************/


  /* CREATE FILENAMES */
  sprintf(mfile,"mdot_%s.dat",root);  
  
  /* INITIALISE PROBLEM */
  /* OPEN ACCRETION RATE OUTPUT FILE */
  mdot=fopen(mfile,"w");
  /* INITIALISE FILE-NAME COUNTER */
  dump_counter = 0;

  /* MAKE GRID */
  /* STAGGERED GRID - Rz & X ARE CELL CENTRES, Rf & XX ARE CELL FACES      */
  /* THE OFFSET (Rf/XX) GRID IS NOT STRICTLY NECESSARY HERE, AS THE FIRST- */
  /* ORDER INTEGRATION SCHEME ONLY REQUIRES THE Rz/X GRID.                 */
  /* IN THIS VERSION THE STAGGERED GRID IS ONLY USED FOR TRACKING MASS     */
  /* CONSERVATION                                                          */

  /* GRID SCALE  - dX=2dr^1/2  */
  /* INNER CELL NOT NECESSARILY AT R=0 */
  /* in & out ARE THE INDICES OF THE BOUNDARY CELLS */
  dX=0.01;
  in=10;
  out=500;
  
  /* MAKE ZONE CENTRED GRID */
  for (i=0; i<out+1; i++)
    {
      X[i]=(dX*((double)i+0.5));   /* UNITS ARE AU^0.5 */
      Rz[i]=(X[i]*X[i])/4.0;
    }
  /* MAKE FACE CENTRED GRID (SOMETIMES CALLED THE "HALF-CELL" GRID) */
  for (i=0; i<out+1; i++)
    {
      XX[i]=(dX*((double)i));  /* AU^0.5 */
      Rf[i]=(XX[i]*XX[i])/4.0;
    }
  
  
  /* PRINT GRID RANGES */
  printf("Zero-torque inner boundary (first non-zero cell) at R=%g\n",Rf[in+1]);
  if (BC_flag==0)
    printf("Zero-torque outer boundary (first non-zero cell) at R=%g\n",Rf[out]);
  if (BC_flag==1)
    printf("Reflecting outer boundary at R=%g\n",Rf[out]);
  
  
  /* SET INITIAL CONDITIONS */
  /* SURFACE DENSITY */
  dR = 0.1;
  for (i=in+1; i<out+1; i++)
    {
      /* SPIKE OF WIDTH dR AT R=1 */
      if ((Rz[i] >= 1.0-dR/2.0) && (Rz[i]< 1.0+dR/2.0))
	sigma[i][0] = 1.0;

      /* NUMERICAL FLOOR FUNCTION */
      /* IMPOSED TO PREVENT CODE CRASH RESULTING FROM sigma <= 0 */
      if (sigma[i][0] < epsilon)
	sigma[i][0] = epsilon;
      
    }
  
  
  /* SET VISCOSITY LAW */
  nu_s = 0.01;
  for (i=in; i<out+1; i++)
    {
      /* CONSTANT VISCOSITY */
      nu[i]= nu_s;  /* AU^2/yr */      
    }
  

  /* SET Y & S VARIABLES */
  for (i=in+1; i<out; i++)
    {
      Y[i][0]=3.0*nu[i]*sqrt(Rz[i])*sigma[i][0];
      S[i][0]=pow(Rz[i],1.5)*sigma[i][0];
    }
  
  
  /* INNER BC */
  sigma[in][0]=0.0;
  S[in][0]=0.0;
  Y[in][0]=0.0;
  
  /* OUTER BC */
  /* ZERO-TORQUE - Sigma=0 */
  if (BC_flag==0)
    {
      sigma[out][0]=0.0;
      S[out][0]=0.0;
      Y[out][0]=0.0;
    }
  /* REFLECTIVE - v_r=0 */
  if (BC_flag==1)
    {
      Y[out][0]=Y[out-1][0];
    }
    
  

  /* COMPUTE ACCRETION RATE */
  M_dot[0]=2.0*(double)M_PI*((Y[in+1][0]-Y[in][0])/dX);
  /* INITIAL MASSES */
  disc_mass=0.0;
  /* ADD UP DISC MASS */
  for (i=1; i<out; i++)
    {
      disc_mass = disc_mass + (2.0*(double)M_PI*Rz[i]*(Rf[i+1]-Rf[i])*sigma[i][0]);
    }
  printf("Initial gas mass = %g\n",disc_mass);
  printf("Initial accretion rate = %g\n",M_dot[0]);

  /* OUTPUT ARRAYS */
  fprintf(mdot,"0.0 %g %g\n",M_dot[0],disc_mass);
  sprintf(sigfile,"sigma_%s_000%i.dat",root,dump_counter);
  sigmaout = fopen(sigfile,"w");
  fprintf(sigmaout,"0.0\n");
  for (i=in; i<out+1; i++)
    {
      fprintf(sigmaout,"%g %g\n",Rz[i],sigma[i][0]);
    }      
  fclose(sigmaout);
  dump_counter++;

 
  
  /* TIME EVOLUTION STARTS HERE */
  /* SET TIME VARIABLES */
  t=1;
  time=0.0;
  time_to_dump=step;
  
  printf("Starting time evolution...\n");

  /* LOOP STARTS HERE */
  while(time <= t_max)
    {

      /* TIME INTEGRATION LOOP */
      /* [0] <-> t ; [1] <-> t+1 */
      /* OVER-WRITES [1] TO [0] AT THE END OF THE LOOP, THEN ITERATES */
            
      /* ------------ STEP 1 - FIND NEW TIMESTEP ------------ */	  
      minstep=1E10;
      /* LOOP GRID */
      for (i=in+1; i<out; i++)
	{
	  /* FROM VISCOSITY */
	  gstep=Courant*dX*dX*(S[i][0]/Y[i][0]);
	  /* FIND MINIMUM */
	  if (gstep < minstep)
	    minstep=gstep;	   
	}

      /* SET TIMESTEP */
      dt = minstep;

      /* UPDATE TIME COUNTER */
      time=time+dt;
      time_to_dump = time_to_dump - dt;

     
      /* ------------ STEP 2 - SURFACE DENSITY ------------ */
      /* LOOP THROUGH GRID AND INTEGRATE */
      for (i=in+1; i<out; i++)
	S[i][1] = S[i][0] + ( ((1.0/(dX*dX)) * (Y[i-1][0]+Y[i+1][0] - 2.0*Y[i][0])) * dt);
	  
      
      /* LOOP GRID AGAIN TO SET OTHER VARIABLES */
      for (i=in+1; i<out; i++)
	{
      	  /* SURFACE DENSITY */
	  sigma[i][1]=S[i][1]/(pow(Rz[i],1.5));
	  /* Y */
	  Y[i][1]=3.0*nu[i]*sqrt(Rz[i])*sigma[i][1];
	}
      

      
      /* ------------ STEP 3 - BOUNDARY CONDITIONS ------------ */	  
      /* SET INNER ZONE TO ZERO: ZERO-TORQUE BC */
      sigma[in][1]=0.0;
      Y[in][1]=0.0;
      S[in][1]=0.0;
      
      /* OUTER BC */
      /* ZERO-TORQUE - Sigma=0 */
      if (BC_flag==0)
	{
	  sigma[out][1]=0.0;
	  S[out][1]=0.0;
	  Y[out][1]=0.0;
	}
      /* REFLECTIVE - v_r=0 */
      if (BC_flag==1)
	{
	  Y[out][1]=Y[out-1][1];
	}
      
      /* ACCRETION RATE AT ORIGIN */
      /* FOR CONSERVATION TRACKING */
      M_dot[1]=2.0*(double)M_PI*((Y[in+1][0]-Y[in][0])/dX);
      /* TRACK ACCRETED MASS (THROUGH INNER BOUNDARY) */
      M_acc = M_acc + M_dot[1]*dt;
      /* TRACK MASS FLOW THROUGH OUTER BOUNDARY */
      M_out = M_out + (dt*2.0*(double)M_PI*((Y[out-1][0]-Y[out][0])/dX));
      /* INTEGRTATE TO FIND TOTAL DISC MASS */
      disc_mass=0.0;
      for (i=in; i<out; i++)
	disc_mass = disc_mass + (2.0*(double)M_PI*Rz[i]*(Rf[i+1]-Rf[i])*sigma[i][1]);



      /* ------------ STEP 4 - DATA I/O & BOOK-KEEPING ------------ */	        
      /* OUTPUT DUMP */
      if (time_to_dump <= 0.0)
	{
	  /* ECHO TO SCREEN */
	  printf("%i t=%g  Mdot=%g  Md=%g\n",t,time,M_dot[1],disc_mass);
	  /* DUMP TO FILES */ 

	  /* CREATE FILENAME */
	  if (dump_counter < 10)
	    sprintf(sigfile,"sigma_%s_000%i.dat",root,dump_counter);
	  if ((dump_counter >= 10) && (dump_counter < 100))
	    sprintf(sigfile,"sigma_%s_00%i.dat",root,dump_counter);
	  if ((dump_counter >= 100) && (dump_counter < 1000))
	    sprintf(sigfile,"sigma_%s_0%i.dat",root,dump_counter);
	  if (dump_counter >=1000)
	    sprintf(sigfile,"sigma_%s_%i.dat",root,dump_counter);

	  /* OPEN SURFACE DENSITY OUTPUT FILE */
	  sigmaout = fopen(sigfile,"w");
	  /* WRITE TO ACCRETION RATE OUTPUT FILE */
	  fprintf(mdot,"%g %g %g\n",time,M_dot[1],disc_mass);
	  /* WRITE TIME TO SURFACE DENSITY FILE */
	  fprintf(sigmaout,"%g\n",time);
	  /* DUMP SURFACE DENSITY */
	  for (i=in; i<out+1; i++)
	    fprintf(sigmaout,"%g  %g\n",Rz[i],sigma[i][1]);
	  /* CLOSE OUTPUT FILE */
	  fclose(sigmaout);

	  /* INCREMENT FILE COUNTER */
	  dump_counter++;
	  /* RESET TIME TO NEXT OUTPUT DUMP */
	  time_to_dump = step - ( time - (round(time/step)*step) );
	  /* FLUSH OUTPUT BUFFERS */
	  fflush(mdot);
	  fflush(stdout);
	}


      /* COPY ARRAYS TO 0 FOR ITERATION */
      for (i=0; i<out+1; i++)
	{
	  S[i][0]=S[i][1];
	  Y[i][0]=Y[i][1];
	  sigma[i][0]=sigma[i][1];
	}


      /* INCREMENT COUNTER */
      t++;

	  
     
    }

  
  
  /* CLOSE ACCRETION RATE DATA FILE */
  fclose(mdot);

  
  exit(EXIT_SUCCESS);
}