/* Calculates time evolution of initial data (u, u_t)=(A*r^5*Gaussian, B*r^5*Gaussian) for cubic focusing Klein-Gordon equation u_tt-u_rr=-u+u^3/r^2 in radial 3d using centered time and centered space discretization. Saves (A,B) in blowup.dat, dispersion.dat, or indecisive.dat depending on outcome */ /***************************************************** ** Parameters **************************************** ******************************************************/ /* minimum amplitude A */ #define MINA -0.0 /* maximum amplitude A */ #define MAXA 0.3 /* stepsize A */ #define INCA 0.01 /* minimum amplitude B */ #define MINB -5.0 /* maximum amplitude B */ #define MAXB 5.0 /* stepsize B */ #define INCB 0.01 /* Grid mesh dr in multiples of DRBASE, see below */ #define DRMULT 30 /* Number of spatial gridpoints */ #define KD 6000 /* Courant constant dt/dr */ #define R 0.9 /* Energy of Q */ #define Qenergy 1.5038 /*********************************************************/ #include #include #include /* define resolution of soliton file */ #define DRBASE 1.0e-4 /* define grid mesh */ #define DR ((double)DRMULT*DRBASE) /* global variable for Q */ double *Q; /****************** initial data *****************/ /* initial position */ double f(double r, double A) { if (r>.5 && r<1) { return A; } else if (r<=0.5) { return 8*A*r*r*r*(1+(0.5-r)*(0.5-r)*(0.5-r)); } else if (r>=1) { return A*(1+(r-1)*(r-1)*(r-1))*exp(-r*r+1); } } /* initial velocity */ double g(double r, double B) { return 10*B*r*r*r*exp(-(r-.5)*(r-.5)); } /***********************************************/ /* import routine for soliton (or other functions) data are assumed to be given with dr=DRBASE caller has to provide pointer to array */ int import(char *file, double *data, long size) { register long k, j; /* loop variables */ double dummy; /* auxiliary variable to store skipped entries */ FILE *fp; /* file pointer */ /* try to open file for reading */ if((fp=fopen(file, "r"))==NULL) return 0; /* no success, return error */ /* loop for reading data */ for(k=0; k1 */ for(j=1; j