text:
/*! \file profile_coord_map.h
    \brief Grid initialization, finding flux coordinates, and getting profile data

    Called by all the top-level functions.
 */

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <Python.h>
#include <gsl/gsl_vector.h>	// needed for 2-d root-finding
#include <gsl/gsl_multiroots.h>
#include <gsl/gsl_linalg.h>
#include "eq_quantity.h"

#ifndef REAL_TYPEDEF
#define REAL_TYPEDEF
#ifndef single			// compiler option determines variable size
typedef double REAL;
#else
typedef float REAL;
#endif
#endif

#ifndef DEBUG			// default debug level: no output to stderr
#define DEBUG 0
#endif

#ifndef DEVELOPMENT		// default development level: use the best of what has been tested to work
#define DEVELOPMENT 0
#endif

#ifndef VERBOSE
#define VERBOSE 0
#endif

#define MAX_N_TRIES 4		// number of times to try to get the root-finder unstuck

// used by desiredrz: constants to choose the desired grid in R,Z space
// used by desiredxyz: constants for desired 3d grid in x,y,z
extern double Xmin,Xmax,Ymin,Ymax,Zmin,Zmax;

// used by getAllProfiles
#define NORM_TO_KEV 4.7894e4      // conversion constant from normalized units to keV (called Tem_00 in fortran files)  Normalization in GTS: Length 100cm, Time Deuterium Cyclotron frequency in 1T B field -- OMEGA_D^-1,
// T normalized st. vth = sqrt(T).  Temperature has the energy dimention.

#define INVCM3_TO_INVM3 1e6	// convert cm^{-3} to m^{-3}

// used by decayNToutsideLCFS and getBoundaryPoints
#define ABOUNDARY 0.9
#define AINSIDE 0.88
#define AZERO 1.2
//#define DECAY_SCALING 0.1*(ABOUNDARY-AINSIDE) // the 0.1 makes it drop off faster

#define MINIMUM(x,y) ((x) < (y) ? (x) : (y))

extern char* NTFileName;

// ---------- local function declarations
int cylToLargeAR(REAL R[],REAL Z[], REAL r[], REAL theta[],int n,REAL R0,REAL Z0);
int largeARtoCyl(REAL r[],REAL theta[], REAL R[], REAL Z[],int n,REAL R0,REAL Z0);
int desiredrz(REAL R[],REAL Z[],REAL R1d[],REAL Z1d[],int nr,int nz,REAL R0,REAL Z0);
int desiredxyz(int nx,int ny,int nz,REAL x1d[],REAL y1d[],REAL z1d[],REAL x[],REAL y[],REAL z[],REAL R0,REAL Z0);
int mesh2dGrid(REAL x2d[],REAL y2d[],REAL x1d[],REAL y1d[],int nx,int ny);
int mesh3dGrid(REAL x3d[],REAL y3d[],REAL z3d[],REAL x1d[],REAL y1d[],REAL z1d[],int nx,int ny,int nz);

int guessatheta(REAL a[],REAL theta[], REAL R[], REAL Z[],int n,REAL R0,REAL Z0);
REAL error_amt(REAL R,REAL Z,REAL a,REAL theta);
void print_state(size_t iter, gsl_multiroot_fsolver * s);
int findcorrectatheta(REAL a[],REAL theta[],REAL Ract[],REAL Zact[],REAL Bm[],REAL Rwant[], REAL Zwant[],int n,int flag[],int location[],int mismatch[]);
int coord_f(const gsl_vector *x,void *params,gsl_vector *f);
int getBoundaryPoints(double** R_bdy,double** Z_bdy,int n_bdy);
void free_BoundaryPoints(double* R_bdy,double* Z_bdy);

int get_mag_axis(REAL coords[]);
int decayNToutsideLCFS(int npts,REAL a[],REAL ne[],REAL Te[],REAL Ti[],int* flag);

int getFluxCoords(int npts,REAL a[],REAL theta[],REAL Bm[],REAL Ract[],REAL Zact[],REAL Rinitial[],REAL Zinitial[],REAL Rwant[],REAL Zwant[],REAL mag_axis_coords[],int flag[],int location[], int mismatch[]);

int getAllProfiles(int npts,REAL Bpol[],REAL BR[],REAL BZ[],REAL Ti[],REAL Te[],REAL P[],REAL ne[],REAL qprofile[],REAL a[],REAL theta[],REAL Rwant[],REAL Zwant[],REAL R_bdy[],REAL Z_bdy[],int InOutFlag[],int location[]);

int getProfile(REAL a,REAL* ne,REAL* Ti,REAL* Te,REAL* a_p,REAL* n_p,REAL* Ti_p,REAL* Te_p,int n_prof );

int adiabaticElectronResponse(int npts,int ntimesteps,REAL ne_tilde[],REAL ne0[],REAL phi[],REAL Te[],int flag[]);

//modified by Lei Shi adding function to find the outside points

//declaration
void inner_check(double R[],double Z[],double R_bdy[],double Z_bdy[],int flag[],int n,int n_bdy,double tol);
//R,Z contains grid point coordinates, R_bdy Z_bdy contains boundary points, flag show inside(with value 1) and outside(with value 0), n the number of grid points, n_bdy the number of boundary points.

//if (R[i],Z[i]) is inside or on the boundary of the contour specified by R_bdy and Z_bdy, flag[i] is set to 1, otherwise 0;

//definition

int horipos_check(double Rc,double Zc, double R1,double Z1,double R2,double Z2,double tol){
  double R=R1+ (Zc-Z1)/(Z2-Z1)*(R2-R1);

  if(fabs(R-Rc)<=tol)
    return 2;
  else if(R<Rc)
    return 1;
  else
    return -1;
}

void inner_check(double R[],double Z[],double R_bdy[],double Z_bdy[],int flag[],int n,int n_bdy,double tol){

  int i;

  for(i=0;i<n;i++){
    double Rc=R[i];
    double Zc=Z[i];

    flag[i]=0;

    int vertpos=0;
    int horipos1=0;
    int horipos2=0;
    int j;
    for(j=0;j<n_bdy;j++){
      if(Z_bdy[j]==Zc){
	if (horipos1==0)
	  {
	    if(fabs(R_bdy[j]-Rc) < tol)
	      horipos1=2;
	    else if(R_bdy[j]<Rc)
	      horipos1=1;
	    else
	      horipos1=-1;
	    continue;
	  }
	else{
	  {
	    if(fabs(R_bdy[j]-Rc) < tol)
	      horipos2=2;
	    else if(R_bdy[j]<Rc)
	      horipos2=1;
	    else
	      horipos2=-1;
	    flag[i]= (horipos1 != horipos2)?1:0;
	    break;
	  }
	}
      }
      else if(Z_bdy[j]<Zc){
	if(vertpos==0){
	  vertpos = -1;
	  continue;
	}
	else if(vertpos==-1)
	  continue;
	else{
	  if(horipos1 == 0){
	    vertpos = -1;
	    horipos1 = horipos_check(Rc,Zc,R_bdy[j],Z_bdy[j],R_bdy[j-1],Z_bdy[j-1],tol);
	    continue;
	  }
	  else{
	    horipos2=horipos_check(Rc,Zc,R_bdy[j],Z_bdy[j],R_bdy[j-1],Z_bdy[j-1],tol);
	    flag[i]= (horipos1 != horipos2)?1:0;
	    break;
	  }

	}

      }
      else{
	if(vertpos==0){
	  vertpos = 1;
	  continue;
	}
	else if(vertpos == 1)
	  continue;
	else{
	  if(horipos1 == 0){
	    vertpos = 1;
	    horipos1 = horipos_check(Rc,Zc,R_bdy[j],Z_bdy[j],R_bdy[j-1],Z_bdy[j-1],tol);
	    continue;
	  }
	  else{
	    horipos2=horipos_check(Rc,Zc,R_bdy[j],Z_bdy[j],R_bdy[j-1],Z_bdy[j-1],tol);
	    flag[i]= (horipos1 != horipos2)?1:0;
	    break;
	  }
	}
      }
    }
  }
}

//end modification

// ---------- data structures
struct coordparams{ 		// structure for 2D root finding
  REAL Rwant;
  REAL Zwant;
};

// ---------- function definitions

//! Convert cylindrical coords (R,Z) to a large aspect ratio approximation of (r,theta)
/*! For n points, where (R0,Z0) is the position of the magnetic axis */
int cylToLargeAR(REAL R[],REAL Z[], REAL r[], REAL theta[],int n,REAL R0,REAL Z0){
  REAL x,z;
  int i;

  for(i=0;i<n;i++){		// here I use: R-R0 = a*cos(theta)
    x = R[i]-R0;		//    and Z = a*sin(theta)
    z = Z[i]-Z0;
    theta[i]=-atan2(z,x);	// poloidal angle (correct quadrant)
    //    r[i] = sqrt(x*x+z*z);	// minor rad.
    r[i] = hypot(x,z);
  }

  return 0;
}

//! Convert a large aspect ratio (r,theta) to cylindrical coords (R,Z)
/*! For n points, where (R0,Z0) is the position of the magnetic axis */
int largeARtoCyl(REAL r[],REAL theta[], REAL R[], REAL Z[],int n,REAL R0,REAL Z0){
  int i;
  for(i=0;i<n;i++){		// here I use: R-R0 = a*cos(theta)
    R[i] = R0 + r[i]*cos(theta[i]);
    Z[i] = Z0 - r[i]*sin(theta[i]);
  }

  return 0;
}

//! Generate a regular, 2D array meshes in (R,Z), and also 1D arrays for each dimension
/*! 2D arrays vary fastest along R (columns) and slowly along Z (rows).
  Also generates 1D arrays R1d, Z1d.
  nr: points in R, nz: # points in z; total # points = nr*nz
  R0,Z0 are the coordinates for the magnetic axis, which is not currently used.
*/
int desiredrz(REAL X[],REAL Y[],REAL X1d[],REAL Y1d[],int nx,int ny,REAL X0,REAL Y0){
  int i;
  REAL dx = (Xmax-Xmin)/(nx-1);
  REAL dy = (Ymax-Ymin)/(ny-1);
  for(i=0;i<ny;i++) Y1d[i] = (REAL)Ymin + i*dy;
  for(i=0;i<nx;i++) X1d[i] = (REAL)Xmin + i*dx;
  mesh2dGrid(X,Y,X1d,Y1d,nx,ny);
  // printf("desiredrz: dr:%g,dz:%g\n",dr,dz);
  // for(i=0;i<nz;i++){
//     Z1d[i] = (REAL)ZMIN + i*dz;
//     for(j=0;j<nr;j++){
//       R[i*nr+j] = (REAL)RMIN + j*dr;
//       Z[i*nr+j] = (REAL)ZMIN + i*dz;
//     }
//   }
  return 0;
}

//! Generate a regular, 3D array meshes in (x,y,z), and also 1D arrays for each dimension
/*! 3D arrays vary fastest along x (columns), then along y (rows), slowest along z (pages).
  x is locally radially outward, y is vertically upward, z is locally toroidal to make a
  right-hand coordinate system.
  nx: points in x direction, ny: points in y direction, nz: points in z direction;
  total # points = nx*ny*nz.
  R0,Z0 are the coordinates for the magnetic axis, which is not currently used.
*/
int desiredxyz(int nx,int ny,int nz,REAL x1d[],REAL y1d[],REAL z1d[],REAL x[],REAL y[],REAL z[],REAL X0,REAL Y0){
  int i,j,k;
  REAL dx = (Xmax-Xmin)/(nx-1);
  REAL dy = (Ymax-Ymin)/(ny-1);
  REAL dz = (Zmax-Zmin)/(nz-1);
#if DEBUG > 1
  fprintf(stderr,"desiredrz: dx:%g,dy:%g,dz:%g\n",dx,dy,dz);
#endif
  for(i=0;i<nz;i++) z1d[i] = (REAL)Zmin + i*dz;
  for(i=0;i<ny;i++) y1d[i] = (REAL)Ymin + i*dy;
  for(i=0;i<nx;i++) x1d[i] = (REAL)Xmin + i*dx;
  mesh3dGrid(x,y,z,x1d,y1d,z1d,nx,ny,nz);
  // for(i=0;i<nz;i++){
//     z1d[i] = (REAL)zMIN + i*dz;
//     for(j=0;j<ny;j++){
//       for(k=0;k<nx;k++){
// 	x[i*ny*nx+j*nx+k] = (REAL)xMIN + k*dx;
// 	y[i*ny*nx+j*nx+k] = (REAL)yMIN + j*dy;
// 	z[i*ny*nx+j*nx+k] = (REAL)zMIN + i*dz;
//       }
//     }
//   }
  return 0;
}

//! Generate 2d mesh grids from 1d grids in each dimension
/*! Takes 1d arrays of coordinates x1d,y1d and generates 2d arrays
    x2d,y2d of nx columns and ny rows, corresponding to the x,y
    coordinate of each point.
*/
int mesh2dGrid(REAL x2d[],REAL y2d[],REAL x1d[],REAL y1d[],int nx,int ny){
  int i,j;
  for(i=0;i<ny;i++){
    for(j=0;j<nx;j++){
      x2d[i*nx+j] = x1d[j];
      y2d[i*nx+j] = y1d[i];
    }
  }

  return 0;
}

//! Generate 3d mesh grids from 1d grids in each dimension
/*! Takes 1d arrays of coordinates x1d,y1d,z1d and generates 3d arrays
    x3d,y3d,z3d of nx columns, ny rows, and nz pages corresponding to the x,y,z
    coordinates of each point.
*/
int mesh3dGrid(REAL x3d[],REAL y3d[],REAL z3d[],REAL x1d[],REAL y1d[],REAL z1d[],int nx,int ny,int nz){
  int i,j,k;
 for(i=0;i<nz;i++){
     for(j=0;j<ny;j++){
      for(k=0;k<nx;k++){
	x3d[i*ny*nx+j*nx+k] = x1d[k];
	y3d[i*ny*nx+j*nx+k] = y1d[j];
	z3d[i*ny*nx+j*nx+k] = z1d[i];
      }
    }
  }

    return 0;
}

//! Given cylindrical coords (R,Z) calculates the initial guess for the flux coords (a,theta)
// first guess as to what the correct a,theta values are for the desired R,Z points
// this is the "large aspect ratio, circ. xsection where I take a = minor radius and
// theta is simply poloidal angle
// then I modify the angle by doing an interpolation
// n is total # of points in 1d array
int guessatheta(REAL a[],REAL theta[], REAL R[], REAL Z[],int n,REAL R0,REAL Z0){
  int i,k;
  int n1dgrid = (int)sqrt(n);
  int ngrid = n1dgrid*n1dgrid;//+10;

#if DEVELOPMENT >= 0

  REAL rwant[n],polwant[n];
  k = cylToLargeAR(R,Z,rwant,polwant,n,R0,Z0);// the minor radius and polodial angles of the R,Z we want

  // --- current method: 2D interpolation for 2 variables:
  //       (rgrid,polgrid) -> a
  //       (rgrid,polgrid) -> theta
  REAL agrid[ngrid],thetagrid[ngrid],Rgrid[ngrid],Zgrid[ngrid];
  REAL rgrid[ngrid],polgrid[ngrid],Bgrid[ngrid];

  // set up a grid to interpolate on
  // use to get an even grid in a,theta
  double da = 1.5/(n1dgrid-1);
  double dtheta = 2.0*M_PI/n1dgrid;
  for(i=0;i<n1dgrid;i++){
    for(k=0;k<n1dgrid;k++){
      agrid[i*n1dgrid+k] = 0.5*da + i*da;
      thetagrid[i*n1dgrid+k] = -M_PI+k*dtheta;
    }
  }

  k = esigetrzb_(Rgrid,Zgrid,Bgrid,agrid,thetagrid,&ngrid); // get Rgrid and Zgrid
  k = cylToLargeAR(Rgrid,Zgrid,rgrid,polgrid,ngrid,R0,Z0);	// get the grid minor radius and poloidal angle

  // interpolate using (rgrid,polgrid) -> (agrid,thetagrid)
  // to get (rwant,thetawant) -> (a,theta)
  k = interp2d(n,ngrid,rwant,polwant,a,theta,rgrid,polgrid,agrid,thetagrid);

#if DEBUG > 0
  fprintf(stderr,"point  agrid   thetagrid  Rgrid  Zgrid  rgrid  polgrid\n");
  for(i=0;i<ngrid;i++)
    fprintf(stderr,"%5d, %g, %g, %g, %g, %g, %g\n",i,agrid[i],thetagrid[i],
	    Rgrid[i],Zgrid[i],rgrid[i],polgrid[i]);
  fprintf(stderr,"point  a   theta  R  Z  rwant  polwant\n");
  for(i=0;i<n;i++){
    fprintf(stderr,"%5d, %g, %g, %g, %g, %g, %g\n",i,a[i],theta[i],
	    R[i],Z[i],rwant[i],polwant[i]);

  }

#endif

    for(i=0;i<n;i++) if(a[i]>2.5) a[i]=2.5;	// limit the maximum radial flux coord

#endif

  // // --- previous method: 1D interpolation of both variables independently
  //   //        polgrid -> (agrid=1,thetagrid)
  //   //        rgrid   -> (agrid,thetagrid=pi/6)

  //   REAL agrid[n1dgrid],thetagrid[n1dgrid],Rgrid[n1dgrid],Zgrid[n1dgrid];
  //   REAL rgrid[n1dgrid],polgrid[n1dgrid],Bgrid[n1dgrid];

  //   // first do the poloidal interpolation
  //   double dtheta = 2*M_PI/n1dgrid;
  //     for(i=0;i<n1dgrid;i++){
  //       agrid[i] = 1;
  //     thetagrid[i] = -M_PI+i*dtheta;
  //     }
  //   k = esigetrzb_(Rgrid,Zgrid,Bgrid,agrid,thetagrid,&n);
  //   k = cylToLargeAR(Rgrid,Zgrid,rgrid,polgrid,n,R0,Z0);	// get the grid minor radius and poloidal angle
  //   k = interpgsl(polwant,theta,polgrid,thetagrid,n,n1dgrid);  // interpolate in poloidal angle:

  //   // then do the radial interpolation
  //   double da = 1.5/(n1dgrid-1);
  //   for(i=0;i<n1dgrid;i++){
  //       agrid[i] = i*da;
  //     thetagrid[i] = M_PI/6;
  //     }
  //   k = esigetrzb_(Rgrid,Zgrid,Bgrid,agrid,thetagrid,&n1dgrid);
  //   k = cylToLargeAR(Rgrid,Zgrid,rgrid,polgrid,n1dgrid,R0,Z0);	// get minor rad. and pol. angle
  //   k = interpgsl(rwant,a,rgrid,agrid,n,n1dgrid); // interpolate radially

  //   // // modification that I discarded
  //   //   for(i=0;i<n;i++){
  //   //     a[i] += 1.5*a[i]*exp(-400.0*theta[i]*theta[i]/(M_PI*M_PI));//a[i]*1.5*fabs(exp(-fabs(theta[i])/M_PI)-exp(-1.0/4.0));//0.5*a[i]*(1-2*fabs(theta[i])/M_PI)*(1-2*fabs(theta[i])/M_PI);
  //   //   }
  //   // --- end method: 1D interpolation of both variables independently

  // --- solution in use previously
#if DEVELOPMENT < 0
  REAL dtheta=2*M_PI/(n-1);// divide flux-coord domain
  REAL x,alarge[n],q[n],Rmap[n],Zmap[n],Bmap[n],qpol[n];//,stheta[n-1];
  // guess a large-aspect ratio, circular cross section
  cylToLargeAR(R,Z,a,theta,n,R0,Z0);
  // for(i=0;i<n;i++){		// here I use: R-R0 = a*cos(theta)
  //     x = R[i]-R0;		//    and Z = a*sin(theta)
  //     theta[i]=-atan2(Z[i],x);	// poloidal angle (correct quadrant)
  //     a[i] = sqrt(x*x+Z[i]*Z[i]);	// minor rad.
  //   }

  // set up interpolation
  // modify large aspect ratio, circ. xsection by interpolating between angles
  // for large minor radius, map theta to actual poloidal angle:
  for(i=0;i<n;i++){//generate an array of evenly-spaced flux coords
    q[i] = -M_PI + dtheta*i;	// to interpolate what qpol (R,Z), q (a,theta)
    alarge[i] = 1;		// is mapped to
  }
  k = esigetrzb_(Rmap,Zmap,Bmap,alarge,q,&n);
  // check to see the poloidal angle these evenly-spaced flux coords mapped to
  for(i=0;i<n;i++){		// here I use: R-R0 = a*cos(theta), Z = a*sin(theta)
    x = Rmap[i]-R0;
    qpol[i]=-atan2(Zmap[i],x);	// poloidal angle (correct quadrant)
    //  rmax[i]=sqrt(x*x+Z[i]*Z[i]);
    //printf("%g,",qpol[i]);
  }
  // now interpolate theta
  REAL thetatmp[n];
  k = interpgsl(theta,thetatmp,qpol,q,n,n);
  for(i=0;i<n;i++) theta[i]=thetatmp[i];      //uncomment to skip interpolation

  // for(i=0;i<n;i++){
  //   a[i] += 1.5*a[i]*exp(-400.0*theta[i]*theta[i]/(M_PI*M_PI));//a[i]*1.5*fabs(exp(-fabs(theta[i])/M_PI)-exp(-1.0/4.0));//0.5*a[i]*(1-2*fabs(theta[i])/M_PI)*(1-2*fabs(theta[i])/M_PI);
  // }  //not worth it; only removed 12 total iterations
  // --- end solution in use previously
#endif

  return 0;

}

//! the 2D function solved by the root finder
int coord_f(const gsl_vector *x,void *params,gsl_vector *f){
  REAL Rwant = ((struct coordparams *) params)->Rwant; // desired R,Z
  REAL Zwant = ((struct coordparams *) params)->Zwant;

  REAL a[1] = {gsl_vector_get (x, 0)}; // current a,theta
  REAL theta[1] = {gsl_vector_get (x, 1)};
  int n=1;
  int k;
  REAL Ract[1],Zact[1],Bm[1]; // actual R,Z at this value of a,theta

  if(!finite(a[0])) a[0] = ABOUNDARY; // in portalr4 this was necessary to eliminate faults
  if(!finite(theta[0])) theta[0] = 0.0;

  k = esigetrzb_(Ract,Zact,Bm,a,theta,&n); // solve for the actual R,Z values

  gsl_vector_set (f, 0, Rwant-Ract[0]); // the error
  gsl_vector_set (f, 1, Zwant-Zact[0]);
  return GSL_SUCCESS;
}

//! returns the distance between desired (R,Z) position and the position mapped to by flux coords (a,theta)
REAL error_amt(REAL R,REAL Z,REAL a,REAL theta){
  REAL Ract,Zact,Bact;
  int single=1;
  esigetrzb_(&Ract,&Zact,&Bact,&a,&theta,&single);

  //  return sqrt((R-Ract)*(R-Ract)+(Z-Zact)*(Z-Zact));
  return hypot(R-Ract,Z-Zact);
}

//! Print current status of the root-finder
void print_state(size_t iter, gsl_multiroot_fsolver * s){
  fprintf (stderr,"iter = %3u x = % .3f % .3f "
	   "f(x) = % .3e % .3e\n",
	   iter,
	   gsl_vector_get (s->x, 0),
	   gsl_vector_get (s->x, 1),
	   gsl_vector_get (s->f, 0),
	   gsl_vector_get (s->f, 1));
}

//! performs 2D root finding to get the (a,theta) for an array of given (R,Z) values
/*! Uses the gsl mulidimensional hybrids root solver ``gsl_multiroot_fsolver_hybrids''.  */
// Modified by Lei Shi
// use function esirz2agq to find (a,theta) for a given (R,Z) inside LCFS, and treat outside points seperately

//write CDF function show the flags

inline void ERR_CI(int e,const char* s)
{
  fprintf(stderr,"error:%s in %s\n",nc_strerror(e),s);
  exit(e);
}

void writeCDF(char* fname,int ntotal, double R[],double Z[],int n_bdy,double R_bdy[],double Z_bdy[],int flag[]){
  int R_dim_id,Z_dim_id,R_bdy_dim_id,Z_bdy_dim_id,flag_dim_id;
  int nr_id,nz_id,ntotal_id,n_bdy_id;
  int R_id,Z_id,R_bdy_id,Z_bdy_id,flag_id;
  int f_id;
  int ecode;
  if(ecode=nc_create(fname,NC_CLOBBER,&f_id))
    ERR_CI(ecode,"create file");
  if(ecode=nc_def_dim(f_id,"r_dim",ntotal,&R_dim_id))
    ERR_CI(ecode,"def r dim");
  if(ecode=nc_def_dim(f_id,"z_dim",ntotal,&Z_dim_id))
    ERR_CI(ecode,"def z dim");
  if(ecode=nc_def_dim(f_id,"r_bdy_dim",n_bdy,&R_bdy_dim_id))
    ERR_CI(ecode,"def r_bdy dim");
  if(ecode=nc_def_dim(f_id,"z__bdy_dim",n_bdy,&Z_bdy_dim_id))
    ERR_CI(ecode,"def z_bdy dim");
  if(ecode=nc_def_dim(f_id,"flag_dim",ntotal,&flag_dim_id))
    ERR_CI(ecode,"def flag dim");
  if(ecode=nc_def_var(f_id,"nr",NC_INT,0,NULL,&nr_id))
    ERR_CI(ecode,"def nr var");
  if(ecode=nc_def_var(f_id,"nz",NC_INT,0,NULL,&nz_id))
    ERR_CI(ecode,"def nz var");
  if(ecode=nc_def_var(f_id,"ntotal",NC_INT,0,NULL,&ntotal_id))
    ERR_CI(ecode,"def ntotal var");
  if(ecode=nc_def_var(f_id,"n_bdy",NC_INT,0,NULL,&n_bdy_id))
    ERR_CI(ecode,"def n_bdy var");
  if(ecode=nc_def_var(f_id,"r",NC_DOUBLE,1,&R_dim_id,&R_id))
    ERR_CI(ecode,"def r var");
  if(ecode=nc_def_var(f_id,"z",NC_DOUBLE,1,&Z_dim_id,&Z_id))
    ERR_CI(ecode,"def z var");
  if(ecode=nc_def_var(f_id,"r_bdy",NC_DOUBLE,1,&R_bdy_dim_id,&R_bdy_id))
    ERR_CI(ecode,"def r_bdy var");
  if(ecode=nc_def_var(f_id,"z_bdy",NC_DOUBLE,1,&Z_bdy_dim_id,&Z_bdy_id))
    ERR_CI(ecode,"def z_bdy var");
  if(ecode=nc_def_var(f_id,"flag",NC_INT,1,&flag_dim_id,&flag_id))
    ERR_CI(ecode,"def flag var");
  if(ecode=nc_enddef(f_id))
    ERR_CI(ecode,"enddef");

  if(ecode=nc_put_var_int(f_id,nr_id,&ntotal))
    ERR_CI(ecode,"put nr var");
  if(ecode=nc_put_var_int(f_id,nz_id,&ntotal))
    ERR_CI(ecode,"put nz var");
  if(ecode=nc_put_var_int(f_id,ntotal_id,&ntotal))
    ERR_CI(ecode,"put ntotal var");
  if(ecode=nc_put_var_int(f_id,n_bdy_id,&n_bdy))
    ERR_CI(ecode,"put n_bdy var");
  if(ecode=nc_put_var_double(f_id,R_id,R))
    ERR_CI(ecode,"put R var");
  if(ecode=nc_put_var_double(f_id,Z_id,Z))
    ERR_CI(ecode,"put Z var");
  if(ecode=nc_put_var_double(f_id,R_bdy_id,R_bdy))
    ERR_CI(ecode,"put R_bdy var");
  if(ecode=nc_put_var_double(f_id,Z_bdy_id,Z_bdy))
    ERR_CI(ecode,"put Z_bdy var");
  if(ecode=nc_put_var_int(f_id,flag_id,flag))
    ERR_CI(ecode,"put flag var");

  if(ecode=nc_close(f_id))
    ERR_CI(ecode,"close file");

  FILE* txtfile;
  txtfile=fopen("./bdy_out.txt","w");
  fprintf(txtfile,"%d\n",n_bdy);
  int i;
  for(i=0;i<n_bdy;i++)
  {
    fprintf(txtfile,"%lf\t%lf\t",R_bdy[i],Z_bdy[i]);
  }
  fclose(txtfile);

}

int findcorrectatheta(REAL a[],REAL theta[],REAL Ract[], REAL Zact[],REAL Bm[], REAL Rwant[], REAL Zwant[],int n,int flag[],int location[], int mismatch[]){
  // find correct a theta coordinates, for inside LCFS points, a theta are obtained by ESIrz2agq function, for outside points, a is assumed to be proportional to the distance from the magnetic axis, and theta is chosen as the same value as the corresponding boundary point. The corresponding boundary point is chosen to be the boundary point closest to the connecting line of magnetic axis and the desired location.
  // Note that in this function, magnetude of B, BR and BZ are also calculated. BR and BZ calculation is not trivial, look into the attached file "Calc_BR_BZ.txt" for more detailed discussion.
  // new version starts

  extern int NBOUNDARY;
  extern double* R_bdy,*Z_bdy;
  fprintf(stderr,"start findcorrecttheta.\n");
  getBoundaryPoints(&R_bdy,&Z_bdy,NBOUNDARY);
  fprintf(stderr,"after getBoundaryPoints.\n");

  double b_axis;
  int single_point=1;
  double tol=1e-9;  //used in inner_check,when given point is near the bdy within tol distance, it's considered on the bdy.
  double dR,dZ;
  double match_tol = 1e-2;

  double theta_guess;
  double theta_bdy;
  double a_guess;
  double a_bdy=ABOUNDARY;
  double zero=0;

  double Z_rel;//the relative Z respect to mag_axis
  double R_rel;//the relative R respect to mag_axis

  double mag_axis_coords[2];
  int loc;
  int sign_ini;
  int sign_cur;

  char CDFfname[20]="./InOutFlags.cdf";
  char logfile[20]="./sl_log.txt";

  FILE* errlog;

  dR = fabs(Rwant[1]-Rwant[0]);
  match_tol *= dR; // absolute tolerance is the percentage tolerance times the finest grid resolution
  fprintf(stderr,"match_tol = %lf\n", match_tol);
  errlog=fopen(logfile,"w+");
  fprintf(errlog,"match_tol = %lf\n", match_tol);
  esigetrzb_(&mag_axis_coords[0],&mag_axis_coords[1],&b_axis,&zero,&zero,&single_point);

  fprintf(stderr,"before inner_check.\n");
  inner_check(Rwant,Zwant,R_bdy,Z_bdy,flag,n,NBOUNDARY,tol);
  fprintf(stderr,"after before write CDF.\n");
  writeCDF(CDFfname,n,Rwant,Zwant,NBOUNDARY,R_bdy,Z_bdy,flag);
  int i;

  /*  for(i=0;i<NBOUNDARY;i++){
    fprintf(errlog,"%lf\t%lf\n",R_bdy[i],Z_bdy[i]);
  }
  fprintf(errlog,"Boundary ends.\n");
  */
  int ierror;
  double b_tmp;

  fprintf(stderr,"Start getting locations.\n");
  for(i=0;i<n;i++){
    if(flag[i]==1){
      //fprintf(stderr,"point %d mapping start,",i);
      int k=ESIrz2agq(&a[i],&theta[i],&Rwant[i],&Zwant[i],&ierror,single_point);
      if(k!=0){
	fprintf(errlog,"Point %d (%lf,%lf) has encountered an error: %d in rz2agq.\n",i,Rwant[i],Zwant[i],ierror);
      }
      esigetrzb_(&Ract[i],&Zact[i],&Bm[i],&a[i],&theta[i],&single_point);
      double err=hypot((Rwant[i]-Ract[i]),(Zwant[i]-Zact[i]));
      double dis=hypot((Rwant[i]-mag_axis_coords[0]),(Zwant[i]-mag_axis_coords[1]));
      // output err and dis for checking
      //fprintf(errlog,"Point %d (%lf,%lf) %lf away from axis has been mapped %lf away\n",i,Rwant[i],Zwant[i],dis,err);
      if( err > match_tol){
	fprintf(errlog,"Point %d (%lf,%lf) %lf away from axis has been mapped %lf away\n",i,Rwant[i],Zwant[i],dis,err);
	mismatch[i] = 1;
	//fprintf(stderr,"Point %d (%lf,%lf) %lf away from axis has been mapped %lf away\n",i,Rwant[i],Zwant[i],dis,err);
	if(fabs(Zwant[i])<0.1)
	{
	  a[i]=0;
	  theta[i]=0;
	}
      }
      location[i]=-1; // inside points do not need location information
    }
    else if(flag[i]==0){
      //      fprintf(stderr,"outter point.\n");
      Bm[i]=b_axis*mag_axis_coords[0]/Rwant[i]; // let B inversely proportional to R, based on B on axis
      Ract[i]=Rwant[i];
      Zact[i]=Zwant[i];
      Z_rel=Zwant[i]-mag_axis_coords[1];
      R_rel=Rwant[i]-mag_axis_coords[0];
      int quarter=0;
      if(R_rel>=0){
	if(Z_rel>=0) quarter=1;
	else quarter =4;
      }
      else{
	if(Z_rel>=0) quarter=2;
	else quarter=3;
      }

      theta_guess=atan((Z_rel/R_rel));
      if(quarter == 2 || quarter == 3)
	theta_guess += M_PI;
      if(quarter == 4)
	theta_guess += 2*M_PI;

      Z_rel=Z_bdy[0]-mag_axis_coords[1];
      R_rel=R_bdy[0]-mag_axis_coords[0];
      theta_bdy=atan((Z_rel/R_rel));

      if(R_rel>=0){
	if(Z_rel>0) quarter=1;
	else quarter =4;
      }
      else{
	if(Z_rel>0) quarter=2;
	else quarter=3;
      }
      if(quarter == 2 || quarter == 3)
	theta_bdy+=M_PI;
      if(quarter == 4)
	theta_bdy += 2*M_PI;

      sign_ini=(theta_guess>=theta_bdy)?1:-1;

      int j;
      for(j=0;j<NBOUNDARY;j++)
	{
	  Z_rel=Z_bdy[j]-mag_axis_coords[1];
	  R_rel=R_bdy[j]-mag_axis_coords[0];
	  theta_bdy=atan((Z_rel/R_rel));

	  if(R_rel>=0){
	    if(Z_rel>0) quarter=1;
	    else quarter =4;
	  }
	  else{
	    if(Z_rel>0) quarter=2;
	    else quarter=3;
	  }
	  if(quarter == 2 || quarter == 3)
	    theta_bdy+=M_PI;
	  if(quarter == 4)
	    theta_bdy += 2*M_PI;//set the theta into (0,2Pi]

	  sign_cur=theta_guess>=theta_bdy?1:-1;
	  if(sign_ini*sign_cur<0){ //if passes the theta we want and is not on the opposite side, then we find the right boundary point
	    location[i]=j;
	    break;
	  }

	}
      if(j==NBOUNDARY)
	location[i]=0;

      Z_rel=Zwant[i]-mag_axis_coords[1];
      R_rel=Rwant[i]-mag_axis_coords[0];
      //      fprintf(stderr,"start guess a\n");
      a_guess=sqrt((Z_rel*Z_rel+R_rel*R_rel)/((R_bdy[loc]-mag_axis_coords[0])*(R_bdy[loc]-mag_axis_coords[0])+(Z_bdy[loc]-mag_axis_coords[1])*(Z_bdy[loc]-mag_axis_coords[1])))*a_bdy;
      a[i]=a_guess;
      theta[i]=2*M_PI/NBOUNDARY*loc;

    }
    else{
      fprintf(stderr,"flag error: point %d at (%lf,%lf) has flag %d.\n",i,Rwant[i],Zwant[i],flag[i]);
      exit(5);
    }
  }
  fprintf(stderr,"finish getting locations.\n");
  fclose(errlog);
  return 0;

  /*old version by Erik
  // declarations for root solver
  const gsl_multiroot_fsolver_type *T;
  gsl_multiroot_fsolver *s;
  int status;			// status of root-finder (notifies if stuck, etc)

  size_t iter = 0;
  int i;
  int ntries;			// num of times for this point that we get stuck
  const size_t ndims = 2;	// 2D coordinate
  struct coordparams p;		// parameters used by ea. instance of solver: Rwant,Zwant
  gsl_multiroot_function f = {&coord_f, ndims, &p}; // root solver setup
  gsl_vector *x = gsl_vector_alloc (ndims); // R,Z values used by root solver

  T = gsl_multiroot_fsolver_hybrids; // solver type
  s = gsl_multiroot_fsolver_alloc (T, 2);

  // loop over all points
  for(i=0;i<n;i++){
    status = GSL_CONTINUE;	// reset root finder status

    // check to see if the nearest point provides a better initial estimate
    if(i>0){
      if(error_amt(Rwant[i],Zwant[i],a[i],theta[i]) > error_amt(Rwant[i],Zwant[i],a[i-1],theta[i-1])){ // if the nearest point is a better guess than the default initial guess, use the previous point
	a[i] = a[i-1];
	theta[i] = theta[i-1];
      }
    }

    // set up this instance of the root finder
    p.Rwant = Rwant[i];		// pass the desired coordinates as parameters
    p.Zwant = Zwant[i];
    gsl_vector_set (x, 0, a[i]);
    gsl_vector_set (x, 1, theta[i]);

    // try a few times with different initial conditions if the solver gets stuck
    for(ntries=0;ntries<MAX_N_TRIES;ntries++){
      if(status == GSL_SUCCESS) break; // in this case, the root-find was successful

      // choose a new initial condition if the solver is stuck
      if(ntries>0){
	#if DEVELOPMENT <= 0
	gsl_vector_set(x,0,0.2);
	gsl_vector_set(x,1,0.25*M_PI+ 0.8*gsl_vector_get(x,1));
	#else
	//	gsl_vector_set(x,0,0.2 + fmod(0.3+gsl_vector_get(x,0),1.2));
	gsl_vector_set(x,0,0.2 + fmod(0.4+gsl_vector_get(x,0),2.0));
	//	gsl_vector_set(x,1,-M_PI + fmod(1.5*M_PI+gsl_vector_get(x,1),2.0*M_PI));
	gsl_vector_set(x,1,-M_PI + fmod(1.25*M_PI+gsl_vector_get(x,1),2.0*M_PI));
	#endif
      }
      gsl_multiroot_fsolver_set (s, &f, x); // initialize the root finder
      iter = 0;			// reset iteration counter
      do{
	iter++;
	status = gsl_multiroot_fsolver_iterate (s);
#if DEBUG > 0
	fprintf(stderr,"point:%d, ntries:%d, status:%d, ",i,ntries,status);
	print_state (iter, s);	// uncomment to give info at each iteration
#endif
	if (status){   // check if solver is stuck
	  break;
	}
	status = gsl_multiroot_test_residual (s->f, 1e-7);
      }while (status == GSL_CONTINUE && iter < 200);

    }

    a[i] = gsl_vector_get(s->x,0);
    theta[i] = gsl_vector_get(s->x,1);
  }

  gsl_multiroot_fsolver_free (s);
  gsl_vector_free (x);

  // make sure all radial flux coords are positive (a>0)
  for(i=0;i<n;i++) a[i]=copysign(a[i],1.0);

  // get the actual R,Z coordinates we've mapped
  status = esigetrzb_(Ract,Zact,Bm,a,theta,&n);

#if VERBOSE > 0
  REAL err_dist;

  //modified by Lei start
  const int NBOUNDARY=4000;
  double R_bdy[NBOUNDARY],Z_bdy[NBOUNDARY];
  getBoundaryPoints(R_bdy,Z_bdy,NBOUNDARY);

  //  for(i=0;i<NBOUNDARY;i++){
  //    fprintf(stdout,"(%lf,%lf),",R_bdy[i],Z_bdy[i]);
  //  }

  for(i=0;i<n;i++){
    err_dist = hypot((Rwant[i]-Ract[i]),(Zwant[i]-Zact[i]));

    if(err_dist > 1e-7)//if the error is unnormally large, try to find a reasonable a and theta for that point, and adjust the B value
      {
	//	fprintf(stderr,"Warning, point %d is mapped %.4g [m] away from requested position.\n",i,err_dist);

	int single_point=1;
	double zero=0;
	double b;         // B field on the axis
	double theta_guess;
	double theta_bdy;
	double a_guess;
	double a_bdy=1;

	double Z_rel;//the relative Z respect to mag_axis
	double R_rel;//the relative R respect to mag_axis

	double mag_axis_coords[2];
	int loc;
	int sign_ini;
	int sign_cur;

	esigetrzb_(&mag_axis_coords[0],&mag_axis_coords[1],&b,&zero,&zero,&single_point);

	Bm[i]=b*mag_axis_coords[0]/Rwant[i]; // let B inversely proportional to R, based on B on axis

	Z_rel=Zwant[i]-mag_axis_coords[1];
	R_rel=Rwant[i]-mag_axis_coords[0];

	theta_guess=atan((Z_rel/R_rel));
	if(theta_guess<=0)
	  theta_guess += 2*M_PI;
	theta_bdy=atan((Z_bdy[0]-mag_axis_coords[1])/(R_bdy[0]-mag_axis_coords[0]));
	if(theta_bdy<=0)
	  theta_bdy+=2*M_PI;
	sign_ini=(theta_guess>=theta_bdy)?1:-1;
	//	fprintf(stderr,"pt[%d]:\n sign_ini = %d\n",i,sign_ini);
	int j;
	for(j=0;j<NBOUNDARY;j++)
	  {
	    theta_bdy=atan((Z_bdy[j]-mag_axis_coords[1])/(R_bdy[j]-mag_axis_coords[0]));
	    if(theta_bdy<=0)
	      theta_bdy += 2*M_PI;//set the theta into [0,2Pi)
	    sign_cur=theta_guess>=theta_bdy?1:-1;
	    if(sign_ini*sign_cur<0 && Zwant[i]*Z_bdy[j]>=0){ //if passes the theta we want and is not on the opposite side, then we find the right boundary point

	      loc=j;
	      break;
	    }

	  }

	a_guess=sqrt((Z_rel*Z_rel+R_rel*R_rel)/((R_bdy[loc]-mag_axis_coords[0])*(R_bdy[loc]-mag_axis_coords[0])+(Z_bdy[loc]-mag_axis_coords[1])*(Z_bdy[loc]-mag_axis_coords[1])))*a_bdy;
	a[i]=a_guess;
	theta[i]=2*M_PI/NBOUNDARY*loc;
	//	if(theta_bdy>M_PI)
	//  fprintf(stderr,"sign_cur=%d\n loc= %d,Rwant=%lf, Zwant=%lf, theta=%lf, Rbd=%lf, Zbd=%lf, thetabd=%lf, a=%lf",sign_cur,loc,Rwant[i],Zwant[i],theta_guess,R_bdy[loc],Z_bdy[loc],theta_bdy,a[i]);
      }
  }

  //Modified by Lei End
#endif

old version ends*/
}

//! Calculates the (R,Z) coordinates of the boundary of the last-closed-flux-surface
/*! Assumes the LCFS is given by the macro ABOUNDARY (set to 1.0 by default) */
int getBoundaryPoints(double** R_bdy,double** Z_bdy,int n_bdy){
  double dtheta = 2*M_PI/(n_bdy-1);
  double a[n_bdy],theta[n_bdy],B_bdy[n_bdy];
  *R_bdy = (double*)xmalloc(sizeof(double)*n_bdy);
  *Z_bdy = (double*)xmalloc(sizeof(double)*n_bdy);
  int i,k;
  for(i=0;i<n_bdy;i++){
    theta[i] = dtheta*i;
    a[i] = ABOUNDARY;
  }

  k=esigetrzb_(*R_bdy,*Z_bdy,B_bdy,a,theta,&n_bdy);

}

void free_BoundaryPoints(double* R_bdy,double* Z_bdy){
  xfree(R_bdy);
  xfree(Z_bdy);
}

//! Find the (R,Z) coordinates of the magnetic axis
/*! Uses the esi functions.
 */
int get_mag_axis(REAL coords[]){
  REAL bm;
  REAL mag_axis = 0;
  int single_point = 1;

  // get the position of the magnetic axis

  //printf("before esigetrzb_\n");
  int k = esigetrzb_(&coords[0],&coords[1],&bm,&mag_axis,&mag_axis,&single_point);
  //printf("after esigetrzb_\n");
  return 0;
}

// modified by Lei Shi to get desired shape outside the LCFS
// The shape function is set to be

// shape1:exponential tan shape
// n(a)=n(1)*exp(-n'(1)/n(1)*tan(PI/2*(a-1)/(a0-1)));

// shape2:cubic shape
// n(a)=c3*a^3+c2*a^2+c1^a+c0

// so that n(1) and n'(1) are unchanged, and n(a0) and n'(a0) are 0;
inline double Shape1(double x, double xb,double x0,double yb,double ypb){
  return yb*exp(ypb/yb*tan(M_PI/2*(x-xb)/(x0-xb)));
}

double Shape2(double x,double xb,double x0,double yb,double ypb){
  double xb2=xb*xb;
  double xb3=xb*xb2;
  double x02=x0*x0;
  double x03=x0*x02;

  double a_data[]=
    {
      1,xb,xb2, xb3,
      1,x0,x02, x03,
      0,1, 2*xb,3*xb2,
      0,1, 2*x0,3*x02
    };

  double b_data[]={yb,0,ypb,0};
  gsl_matrix_view m= gsl_matrix_view_array(a_data,4,4);
  gsl_vector_view b= gsl_vector_view_array(b_data,4);
  gsl_vector *c = gsl_vector_alloc(4);
  int s;
  gsl_permutation *p=gsl_permutation_alloc(4);
  gsl_linalg_LU_decomp(&m.matrix,p,&s);
  gsl_linalg_LU_solve(&m.matrix,p,&b.vector,c);

  double C[4];
  int i;
  for(i=0;i<4;i++)
    C[i]=gsl_vector_get(c,i);

  return C[0]+C[1]*x+C[2]*x*x+C[3]*x*x*x;
}

inline double Decay(int shape,double a,double a_bdy,double a0,double y_bdy,double yp_bdy){
  switch(shape){
  case 1:
    return Shape1(a,a_bdy,a0,y_bdy,yp_bdy);
  case 2:
    return Shape2(a,a_bdy,a0,y_bdy,yp_bdy);
  default:
    fprintf(stderr,"Wrong decay shape number, check the notes in profile_coord_map.h before function decayNToutsideLCFS for more info.\n");
  }
}

//! Cause density and temperature to exponentially decay outside the LCFS
/*! Use the value of density, Ti, and Te at the Last-Closed-Flux-Surface
  and force these quantities to exponentially decay outside this surface.
  This is done so there is no abrupt increase in these quantities that might
  frustrate the paraxial solver in the reflectometry code.
*/

int decayNToutsideLCFS(int npts,REAL a[],REAL ne[],REAL Te[],REAL Ti[],int* flag){
  REAL abound = ABOUNDARY; // (a,theta) on the boundary that we use
  REAL ain = AINSIDE;	// (a,theta) inside the LCFS we use
  // (* to find the exponential drop off *) in our case, to find n'(1)

  REAL a0 = AZERO; // the flux coord where we set n & T to be zero

  //load the equilibrium profile
  int** n_prof=(int**)PyMem_Malloc(sizeof(int*));
  REAL **a_p=(REAL**)PyMem_Malloc(sizeof(REAL*));
  REAL **n_p=(REAL**)PyMem_Malloc(sizeof(REAL*));
  REAL **Ti_p=(REAL**)PyMem_Malloc(sizeof(REAL*));
  REAL **Te_p=(REAL**)PyMem_Malloc(sizeof(REAL*));
  int k;
  k=loadNTProfiles(NTFileName,n_prof,a_p,n_p,Ti_p,Te_p);
  if(k!=0){
    fprintf(stderr,"Loading NTProfile error %d.\n",k);
    exit(1);
  }

  // values of ne,te,ti at boundary and just inside LCFS
  REAL nebound,tibound,tebound,nein,tin,tein;
  getProfile(abound,&nebound,&tibound,&tebound,*a_p,*n_p,*Ti_p,*Te_p,**n_prof);
  getProfile(ain,&nein,&tin,&tein,*a_p,*n_p,*Ti_p,*Te_p,**n_prof);

  nebound *=INVCM3_TO_INVM3;
//  tibound *=1;
//  tebound *=1;
  nein *=INVCM3_TO_INVM3;
//  tin *=1;
//  tein *=1;

  // the 1/e distance in radial flux coord
  // the scale factor is to make it drop off faster
  //  REAL lambda_ne = -fabs(DECAY_SCALING*nebound/(nein-nebound));
  //  REAL lambda_ti = -fabs(DECAY_SCALING*tibound/(tin-tibound));
  //  REAL lambda_te = -fabs(DECAY_SCALING*tebound/(tein-tebound));

  // the direvitives are obtained based on linear estimation

  REAL ne_prime = (nebound-nein)/(abound-ain);
  REAL ti_prime = (tibound-tin)/(abound-ain);
  REAL te_prime = (tebound-tein)/(abound-ain);

  fprintf(stderr,"n'=%lf,t'=%lf",ne_prime,te_prime);

  // update the values in the ne,te,ti arrays
//  int i;
//  for(i=0;i<npts;i++){
//    if(a[i]>1.0){
//    ne[i] = nebound*exp((a[i]-abound)/lambda_ne);
//    Te[i] = tebound*exp((a[i]-abound)/lambda_te);
//    Ti[i] = tibound*exp((a[i]-abound)/lambda_ti);
//    }
  int i;
  for(i=0;i<npts;i++){
    if(flag[i] == 0 ){
      if(a[i] <= a0){
	ne[i]=Decay(2,a[i],abound,a0,nebound,ne_prime);
	Te[i]=Decay(2,a[i],abound,a0,tebound,te_prime);
	Ti[i]=Decay(2,a[i],abound,a0,tibound,ti_prime);
	//	fprintf(stderr,"inside a[%d]=%lf, ne=%lf, Te=%lf , Ti=%lf\n",i,a[i],ne[i],Te[i],Ti[i]);
      }
      else{
	ne[i]=Te[i]=Ti[i]=0;  // if the point is outside a0, set all to be 0
	//	fprintf(stderr,"outside a[%d]=%lf, ne=%lf, Te=%lf , Ti=%lf\n",i,a[i],ne[i],Te[i],Ti[i]);
      }
    }
  }
  return 0;
}

//! Runs all the functions to map from cylindrical (R,Z) coords to flux coords
/*! Given the number of points npts, the desired cylindrical coords (Rwant,Zwant)
    and the (R,Z) coordinate array of the magnetic axis position "mag_axis_coords",
    gives the initial guess (Rinitial,Zinitial), the flux coords (a,theta), and
    the actual (Ract,Zact) these flux coords map to.
 */
int getFluxCoords(int npts,REAL a[],REAL theta[],REAL Bm[],REAL Ract[],REAL Zact[],REAL Rinitial[],REAL Zinitial[],REAL Rwant[],REAL Zwant[],REAL mag_axis_coords[],int flag[],int location[], int mismatch[]){
  // field-line coords: a: flux (radial), theta: angle (poloidal), Bm: |B|
  // (Rwant,Zwant): the R,Z coordinatew we want to map to
  // (Rinitial,Zinitial): the R,Z coords of our initial guess
  // (Ract,Zact): actual R,Z coordinates we end up with
  int k;
  //  k = guessatheta(a,theta,Rwant,Zwant,npts,mag_axis_coords[0],mag_axis_coords[1]); // initial guess for the flux coords
  //  k = esigetrzb_(Rinitial,Zinitial,Bm,a,theta,&npts); // map to cylindrical coords

  fprintf(stderr, "before findcorrectatheta.\n");
  k = findcorrectatheta(a,theta,Ract,Zact,Bm,Rwant,Zwant,npts,flag,location,mismatch); // use esi provided interpolation to get flux coords

  return k;
}

//! Gets all the profile data at the specified flux coordinates
/*! Given the number of points npts, the flux coords (a,theta), gets the
    values of: Bpol, Ti, Te, P, ne, and the q profile at these points.
 */
extern char* NTFileName;

int getProfile(REAL a,REAL* ne,REAL* Ti,REAL* Te,REAL* a_p,REAL* n_p,REAL* Ti_p,REAL* Te_p,int n_prof){ //get ne,Ti,Te for a specified flux surface psi=a, linear interpolation based on the loaded ne,Ti,Te profile

  interpxy(&a,ne,a_p,n_p,1,n_prof);
  interpxy(&a,Ti,a_p,Ti_p,1,n_prof);
  interpxy(&a,Te,a_p,Te_p,1,n_prof);
  return 0;
}

int getAllProfiles(int npts,REAL Bpol[], REAL BR[],REAL BZ[] ,REAL Ti[],REAL Te[],REAL P[],REAL ne[],REAL qprofile[],REAL a[],REAL theta[], REAL Rwant[], REAL Zwant[],REAL R_bdy[],REAL Z_bdy[],int InOutFlag[], int location[]){
  // used by esilink2c and esiget2dfunctions
  REAL *F = (REAL *)PyMem_Malloc(npts*sizeof(REAL));
  REAL *Fa = (REAL *)PyMem_Malloc(npts*sizeof(REAL));
  REAL *gFa = (REAL *)PyMem_Malloc(npts*sizeof(REAL));
  REAL *gFaa = (REAL *)PyMem_Malloc(npts*sizeof(REAL));
  REAL *gYa = (REAL *)PyMem_Malloc(npts*sizeof(REAL));
  REAL *gYaa = (REAL *)PyMem_Malloc(npts*sizeof(REAL));
  REAL *Ta = (REAL *)PyMem_Malloc(npts*sizeof(REAL));
  REAL *Pa = (REAL *)PyMem_Malloc(npts*sizeof(REAL));
  REAL *r = (REAL *)PyMem_Malloc(npts*sizeof(REAL));
  REAL *ra = (REAL *)PyMem_Malloc(npts*sizeof(REAL));
  REAL *rq = (REAL *)PyMem_Malloc(npts*sizeof(REAL));
  REAL *z = (REAL *)PyMem_Malloc(npts*sizeof(REAL));
  REAL *za = (REAL *)PyMem_Malloc(npts*sizeof(REAL));
  REAL *zq = (REAL *)PyMem_Malloc(npts*sizeof(REAL));
  REAL *B = (REAL *)PyMem_Malloc(npts*sizeof(REAL));
  REAL *Ba = (REAL *)PyMem_Malloc(npts*sizeof(REAL));
  REAL *Bq = (REAL *)PyMem_Malloc(npts*sizeof(REAL));
  REAL *gh = (REAL *)PyMem_Malloc(npts*sizeof(REAL));
  REAL *gha = (REAL *)PyMem_Malloc(npts*sizeof(REAL));
  REAL *ghq = (REAL *)PyMem_Malloc(npts*sizeof(REAL));

  int k = esilink2c_(F,Fa,gFa,gFaa,gYa,gYaa,Ti,Ta,P,Pa,r,ra,rq,z,za,zq,B,Ba,Bq,gh,gha,ghq);
  k = esiget2dfunctions_(a,theta,&npts);

  // get q profile
  int i;
  for(i=0;i<npts;i++){// reverse of method documented in escZ_march03.c
    if(a[i]==0.0){		// this method is used in setup_v2.F90
      qprofile[i] = -gFaa[i]/gYaa[i]; // the magnetic axis is a special case
    }else qprofile[i]=-gFa[i]/gYa[i];
  }

  int** n_prof=(int**)PyMem_Malloc(sizeof(int*));
  REAL **a_p=(REAL**)PyMem_Malloc(sizeof(REAL*));
  REAL **n_p=(REAL**)PyMem_Malloc(sizeof(REAL*));
  REAL **Ti_p=(REAL**)PyMem_Malloc(sizeof(REAL*));
  REAL **Te_p=(REAL**)PyMem_Malloc(sizeof(REAL*));

  printf("NTFileName before loadNTProfiles: %s\n",NTFileName);
  k=loadNTProfiles(NTFileName,n_prof,a_p,n_p,Ti_p,Te_p);
  if(k!=0){
    fprintf(stderr,"Loading NTProfile error %d.\n",k);
    exit(1);
  }

  // get density and temperature profiles
  for(i=0;i<npts;i++){
    if(InOutFlag[i]==1){
      getProfile(a[i],&ne[i],&Ti[i],&Te[i],*a_p,*n_p,*Ti_p,*Te_p,**n_prof);
      Ti[i] *= 100;
      Te[i] *= 100;
      ne[i] *= INVCM3_TO_INVM3;
    }
    else{
      ne[i]=0;
      Ti[i]=0;
      Te[i]=0;
    }
  }

  //get BR, BZ for new FWR3D
  // the basic idea is using flux coordinates to express BR = B dot grad(R).
  // look into attached pdf file for more details about calculating these two components.

  fprintf(stderr,"before getting BR,BZ. \n");
  //get magnetic axis coordinates
  double mag_axis_coords[2];
  get_mag_axis(mag_axis_coords);

  for (i=0;i<npts;i++){
    if(InOutFlag[i]==1){

      // for inside points, all relavent quantites can be gotten from ESI equilibrium reading functions defined in eq_quantity.h

      //fprintf(stderr,"inside point [%d].\n",i);
      initialize_esi_out(&a[i],&theta[i],1);
      //fprintf(stderr,"esi_out initialized.\n");
      //fprintf(stderr,"a=%lf,theta=%lf,n=%d.\n",esi_output.aaa[0],esi_output.the[0],esi_output.n);
      load_esi_quantities();
      //fprintf(stderr,"esi_quantities loaded\n");
      double** metric = get_metric_elements();
      //fprintf(stderr,"metric_elements got.\n");
      double gaa = metric[0][0];
      double gtt = metric[1][0];
      double gat = metric[3][0];
      double* J = get_Jacobian();
      //fprintf(stderr,"Jacobian got.\n");
      double dpsi_da = -esi_output.gYaj[0];
      double dr_dth = esi_output.rqj[0];
      double dr_da = esi_output.raj[0];
      double dz_dth = esi_output.zqj[0];
      double dz_da = esi_output.zaj[0];

      // now, define some intermediate quantities:
      double mod_grad_R = sqrt(dr_da*dr_da*gaa + dr_dth*dr_dth*gtt + 2*dr_da*dr_dth*gat);
      double mod_grad_Z = sqrt(dz_da*dz_da*gaa + dz_dth*dz_dth*gtt + 2*dz_da*dz_dth*gat);

      // calculate BR, BZ at this point
      //fprintf(stderr,"before assigning BR,BZ. \n");
      BR[i] = -dpsi_da*dr_dth/(mod_grad_R*J[0]);
      BZ[i] = -dpsi_da*dz_dth/(mod_grad_Z*J[0]);

      // clean up
      //fprintf(stderr,"before cleaning up.");
      free_all(metric,J);
    }
    else{
      //if the point is outside of the LCFS, we estimate the BR, BZ using the information gotten from the corresponding points ON the LCFS. Assume Bpol is decaying proportionally to the distance away from the magnetic axis.
      //fprintf(stderr,"outside point[%d].\n",i);
      extern double* R_bdy,*Z_bdy;
      extern int NBOUNDARY;
      int loc = location[i];
      //fprintf(stderr,"location got.\n");

      double theta_loc = 2*M_PI/(NBOUNDARY-1) * loc; //Note that the convention is that the first and last boundary points corresponds to the same location(theta =0 and 2*pi), so boundary points form a closed circle.
      double a_loc = ABOUNDARY;
      double Rb = R_bdy[loc];
      double Zb = Z_bdy[loc];

      //fprintf(stderr,"Rb,Zb got.\n");

      // for the boundary point, do the same thing to obtain BR, BZ there.
      initialize_esi_out(&a_loc,&theta_loc,1);
      //fprintf(stderr,"esi_out initialized.\n");
      load_esi_quantities();
      //fprintf(stderr,"esi_quantities loaded. \n");
      double** metric = get_metric_elements();
      //fprintf(stderr,"metric elements got.\n");
      double gaa = metric[0][0];
      double gtt = metric[1][0];
      double gat = metric[3][0];
      double* J = get_Jacobian();
      //fprintf(stderr,"Jacobian got.\n");
      double dpsi_da = -esi_output.gYaj[0];
      double dr_dth = esi_output.rqj[0];
      double dr_da = esi_output.raj[0];
      double dz_dth = esi_output.zqj[0];
      double dz_da = esi_output.zaj[0];

      // now, define some intermediate quantities:
      double mod_grad_R = sqrt(dr_da*dr_da*gaa + dr_dth*dr_dth*gtt + 2*dr_da*dr_dth*gat);
      double mod_grad_Z = sqrt(dz_da*dz_da*gaa + dz_dth*dz_dth*gtt + 2*dz_da*dz_dth*gat);

      // calculate BR, BZ at this point
      double BR_bdy = -dpsi_da*dr_dth/(mod_grad_R*J[0]);
      double BZ_bdy = -dpsi_da*dz_dth/(mod_grad_Z*J[0]);
      //fprintf(stderr,"BR,BZ bdy got.\n");
      // clean up
      free_all(metric,J);

      // now calculate the BR, BZ at the outside point
      double R = Rwant[i];
      double Z = Zwant[i];

      //calculate the distances between relavent points
      double dis_want = hypot(R-mag_axis_coords[0],Z-mag_axis_coords[1]);
      double dis_bdy = hypot(Rb-mag_axis_coords[0],Zb-mag_axis_coords[1]);

      //get BR, BZ based on the ratio of distances to magnetic axis
      BR[i] = BR_bdy * dis_bdy / dis_want;
      BZ[i] = BZ_bdy * dis_bdy / dis_want;

    }
  }

  //calculate Bpol
  for(i=0;i<npts;i++){
    Bpol[i] = sqrt(BR[i]*BR[i] + BZ[i]*BZ[i]);
  }

  fprintf(stderr,"after getting BR,BZ. \n");

  // free arrays that are no longer needed:
  clearNTProfiles(n_prof,a_p,n_p,Ti_p,Te_p);
  PyMem_Free(F);
  PyMem_Free(Fa);
  PyMem_Free(gFa);
  PyMem_Free(gFaa);
  PyMem_Free(gYa);
  PyMem_Free(gYaa);
  PyMem_Free(Ta);
  PyMem_Free(Pa);
  PyMem_Free(r);
  PyMem_Free(ra);
  PyMem_Free(rq);
  PyMem_Free(z);
  PyMem_Free(za);
  PyMem_Free(zq);
  PyMem_Free(B);
  PyMem_Free(Ba);
  PyMem_Free(Bq);
  PyMem_Free(gh);
  PyMem_Free(gha);
  PyMem_Free(ghq);

  return 0;
}

//! Determine the density fluctuation assuming adiabatic electron response to potential fluctuations.
/*! Given the number of points npts, the background density (ne0), electron temperature (Te),
    and potential fluctuation (phi), computes the density fluctuation (ne_tilde) and the total
    density (ne), using: ne_tilde = ne0*phi/Te, and ne = ne0 + ne_tilde
 */
int adiabaticElectronResponse(int npts,int ntimesteps,REAL ne_tilde[],REAL ne0[],REAL phi[],REAL Te[],int flag[]){
    int i,j;
    for(i=0;i<ntimesteps;i++){
      for(j=0;j<npts;j++){
	if(flag[j]==1){
	  double* phi_this = &phi[i*npts+j];
	  *phi_this *= NORM_TO_KEV;
	  ne_tilde[i*npts+j] =  (*phi_this)/Te[j];
	}
	else
	  ne_tilde[i*npts+j] = 0;
      }
    }
    return 0;
  }

max_stars_repo_name: LeiShi/Synthetic-Diagnostics-Platform


text:
/*
 * author: Pierre Schnizer
 * created: March 2004
 * file: pygsl/src/initmodule.c
 * $Id: initmodule.c,v 1.14 2008/10/25 20:45:04 schnizer Exp $
 *
 * Changes:
 *     7. October 2003:
 *     Removed the error handler from this file. It is now in
 *     Lib/error_helpers.c Each module must call init_pygsl()
 *     in its init routine. This is necessary to support platforms
 *     where the gsl library is statically linked to the various
 *     modules.
 */
#define  _PyGSL_API_MODULE 1
#include <pygsl/intern.h>
#include <pygsl/utils.h>
#include <pygsl/error_helpers.h>

#include <gsl/gsl_rng.h>
#include <gsl/gsl_ieee_utils.h>
#include <gsl/gsl_version.h>

/* Taken from Modules/getbuildinfo */
#ifndef DATE
#ifdef __DATE__
#define DATE __DATE__
#else
#define DATE "xx/xx/xx"
#endif
#endif

#ifndef TIME
#ifdef __TIME__
#define TIME __TIME__
#else
#define TIME "xx:xx:xx"
#endif
#endif
/* End */

/*
 * Used as a buffer to generate error messages.
 */
static char pygsl_error_str[512];
#include "profile.c"
#include "error_helpers.c"
#include "general_helpers.c"
#include "complex_helpers.c"
#include "block_helpers.c"
#include "function_helpers.c"
#include "rng_helpers.c"

static PyObject * debuglist = NULL;

static int
PyGSL_register_debug_flag(int * ptr, const char * module_name)
{
#if DEBUG == 1
     PyObject * cobj;
     FUNC_MESS_BEGIN();
     if ((cobj = PyCObject_FromVoidPtr((void *) ptr, NULL)) == NULL){
	  fprintf(stderr, "Could not create PyCObject for ptr %p to debug flag for module %s\n",
		  (void * ) ptr, module_name);
	  return GSL_EFAILED;
     }
     DEBUG_MESS(2, "Registering ptr %p for module %s", (void *) ptr,  module_name);
     if(PyList_Append(debuglist, cobj) != 0){
	  return GSL_EFAILED;
     }
     *ptr = pygsl_debug_level;
     FUNC_MESS_END();
     return GSL_SUCCESS;
#endif
     PyGSL_ERROR("Why does it try to register the debug level callback? Should"
	       " use the compile time value DEBUG!", GSL_ESANITY);
}

static PyObject *
PyGSL_set_debug_level(PyObject *self, PyObject *args)
{

     FUNC_MESS_BEGIN();
#if DEBUG == 1
     PyObject *o;
     int tmp, i, max, *ptr;
     if(!PyArg_ParseTuple(args, "i", &tmp))
	  return NULL;
     if(tmp >= 0 && tmp <PyGSL_DEBUG_MAX)
	  ;
     else
	  PyGSL_ERROR_VAL("Only accept debug levels between 0 and PyGSL_DEBUG_MAX", GSL_EINVAL, NULL);

     pygsl_debug_level = tmp;
     max = PySequence_Size(debuglist);
     DEBUG_MESS(3, "Setting debug level to %d for %d modules", pygsl_debug_level, max);
     for(i = 0; i < max; ++i){
	  if((o = PySequence_GetItem(debuglist, i)) == NULL){
	       fprintf(stderr, "In file %s at line %d; Could not get element %d\n",
		       __FILE__, __LINE__, i);
	       continue;
	  }
	  ptr = (int *)PyCObject_AsVoidPtr(o);
	  DEBUG_MESS(2, "Setting info ptr %p", (void *) ptr);
	  *ptr = tmp;
     }
     Py_INCREF(Py_None);
     FUNC_MESS_END();
     return Py_None;
#else
     PyGSL_ERROR_NULL("PyGSL was not compiled with DEBUG = 1; Can not set DEBUG level!", GSL_EUNIMPL);
#endif
}

static PyObject *
PyGSL_get_debug_level(PyObject *self, PyObject *args)
{
     int tmp;
#if DEBUG == 1
     tmp = (int) pygsl_debug_level;
#else
     tmp = DEBUG;
#endif
     return PyInt_FromLong(tmp);
}

static void * _PyGSL_API[PyGSL_NENTRIES_NUM];

static PyMethodDef initMethods[] = {
     {"get_debug_level", PyGSL_get_debug_level, METH_NOARGS, NULL},
     {"set_debug_level", PyGSL_set_debug_level, METH_VARARGS, NULL},
     {"vector_transform_counter",  PyGSL_get_vector_transform_counter ,METH_NOARGS, NULL},
     {"matrix_transform_counter",  PyGSL_get_matrix_transform_counter ,METH_NOARGS, NULL},
     {"complex_transform_counter", PyGSL_get_complex_transform_counter,METH_NOARGS, NULL},
     {"float_transform_counter",   PyGSL_get_float_transform_counter  ,METH_NOARGS, NULL},
     {"register_exceptions",       PyGSL_register_exceptions,          METH_VARARGS,NULL},
     {"register_warnings",         PyGSL_register_warnings,            METH_VARARGS,NULL},
     {NULL,     NULL, 0, NULL}        /* Sentinel */
};

static void
PyGSL_init_api(void)
{
     int i;

     for(i=0;i<PyGSL_NENTRIES_NUM; ++i){
	  _PyGSL_API[i] = NULL;
     }
     _PyGSL_API[PyGSL_api_version_NUM                          ] = (void *) PyGSL_API_VERSION;
     _PyGSL_API[PyGSL_RNG_ObjectType_NUM                       ] = NULL;
     _PyGSL_API[PyGSL_error_flag_NUM                           ] = (void *) &PyGSL_error_flag;
     _PyGSL_API[PyGSL_error_flag_to_pyint_NUM                  ] = (void *) &PyGSL_error_flag_to_pyint;
     _PyGSL_API[PyGSL_add_traceback_NUM                        ] = (void *) &PyGSL_add_traceback;
     _PyGSL_API[PyGSL_module_error_handler_NUM                 ] = (void *) &PyGSL_module_error_handler;

     _PyGSL_API[PyGSL_error_string_for_callback_NUM            ] = (void *) & PyGSL_set_error_string_for_callback        ;
     _PyGSL_API[PyGSL_pyfloat_to_double_NUM                    ] = (void *) & PyGSL_pyfloat_to_double                    ;
     _PyGSL_API[PyGSL_pylong_to_ulong_NUM                      ] = (void *) & PyGSL_pylong_to_ulong                      ;
     _PyGSL_API[PyGSL_pylong_to_uint_NUM                       ] = (void *) & PyGSL_pylong_to_uint                       ;
     _PyGSL_API[PyGSL_check_python_return_NUM                  ] = (void *) & PyGSL_check_python_return                  ;
     _PyGSL_API[PyGSL_clear_name_NUM                           ] = (void *) & PyGSL_clear_name                           ;

     _PyGSL_API[PyGSL_PyComplex_to_gsl_complex_NUM             ] = (void *) & PyGSL_PyComplex_to_gsl_complex             ;
     _PyGSL_API[PyGSL_PyComplex_to_gsl_complex_float_NUM       ] = (void *) & PyGSL_PyComplex_to_gsl_complex_float       ;
     _PyGSL_API[PyGSL_PyComplex_to_gsl_complex_long_double_NUM ] = (void *) & PyGSL_PyComplex_to_gsl_complex_long_double ;

     _PyGSL_API[PyGSL_stride_recalc_NUM                        ] = (void *) & PyGSL_stride_recalc                        ;
     _PyGSL_API[PyGSL_PyArray_new_NUM                          ] = (void *) & PyGSL_New_Array                            ;
     _PyGSL_API[PyGSL_PyArray_copy_NUM                         ] = (void *) & PyGSL_Copy_Array                           ;
/*
     _PyGSL_API[PyGSL_PyArray_prepare_gsl_vector_view_NUM      ] = (void *) & PyGSL_PyArray_prepare_gsl_vector_view      ;
     _PyGSL_API[PyGSL_PyArray_prepare_gsl_matrix_view_NUM      ] = (void *) & PyGSL_PyArray_prepare_gsl_matrix_view      ;
*/
     _PyGSL_API[PyGSL_PyArray_generate_gsl_vector_view_NUM     ] = (void *) & PyGSL_PyArray_generate_gsl_vector_view     ;
     _PyGSL_API[PyGSL_PyArray_generate_gsl_matrix_view_NUM     ] = (void *) & PyGSL_PyArray_generate_gsl_matrix_view     ;
     _PyGSL_API[PyGSL_copy_pyarray_to_gslvector_NUM            ] = (void *) & PyGSL_copy_pyarray_to_gslvector            ;
     _PyGSL_API[PyGSL_copy_pyarray_to_gslmatrix_NUM            ] = (void *) & PyGSL_copy_pyarray_to_gslmatrix            ;
     _PyGSL_API[PyGSL_copy_gslvector_to_pyarray_NUM            ] = (void *) & PyGSL_copy_gslvector_to_pyarray            ;
     _PyGSL_API[PyGSL_copy_gslmatrix_to_pyarray_NUM            ] = (void *) & PyGSL_copy_gslmatrix_to_pyarray            ;

     _PyGSL_API[PyGSL_gsl_rng_from_pyobject_NUM                ] = (void *) & PyGSL_gsl_rng_from_pyobject                ;
     _PyGSL_API[PyGSL_function_wrap_helper_NUM                 ] = (void *) & PyGSL_function_wrap_helper                 ;
     _PyGSL_API[PyGSL_register_debug_flag_NUM                  ] = (void *) & PyGSL_register_debug_flag                  ;
     _PyGSL_API[PyGSL_vector_or_double_NUM                     ] = (void *) & PyGSL_vector_or_double                     ;
     _PyGSL_API[PyGSL_warning_NUM                              ] = (void *) & PyGSL_warning                              ;
     _PyGSL_API[PyGSL_pyint_to_int_NUM                         ] = (void *) & PyGSL_pyint_to_int                         ;
     _PyGSL_API[PyGSL_vector_check_NUM                         ] = (void *) & PyGSL_vector_check                         ;
     _PyGSL_API[PyGSL_matrix_check_NUM                         ] = (void *) & PyGSL_matrix_check                         ;
     _PyGSL_API[PyGSL_array_check_NUM                          ] = (void *) & PyGSL_array_check                          ;

}

DL_EXPORT(void) initinit(void)
{
  PyObject *m = NULL, *d = NULL, *version=NULL, *date=NULL, *api = NULL;

  m = Py_InitModule("pygsl.init", initMethods);
  import_array();

  if(m == NULL){
       fprintf(stderr, "I could not init pygsl.init!");
       return;
  }

  d = PyModule_GetDict(m);
  if(d == NULL){
       fprintf(stderr, "I could not get the module dict for  pygsl.init!");
       return;
  }

  PyGSL_init_api();

  if(PyGSL_init_errno() != 0){
       PyErr_SetString(PyExc_ImportError, "Failed to init errno handling!");
  }
  PyGSL_API = _PyGSL_API;
  PyGSL_SET_ERROR_HANDLER();

  api = PyCObject_FromVoidPtr((void *) PyGSL_API, NULL);
  assert(api);
  if (PyDict_SetItemString(d, "_PYGSL_API", api) != 0){
       PyErr_SetString(PyExc_ImportError,
		       "I could not add  _PYGSL_API!");
       return;
  }

  version = PyString_FromString(GSL_VERSION);
  if(version == NULL){
       fprintf(stderr, "I could not create the version string for pygsl.init!");
       return;
  }
  if(PyDict_SetItemString(d, "compiled_gsl_version", version) != 0){
       fprintf(stderr, "I could not add the compile version string to the module dict of pygsl.init!");
       return;
  }

  version = PyString_FromString(gsl_version);
  if(version == NULL){
       fprintf(stderr, "I could not create the version string for pygsl.init!");
       return;
  }
  if(PyDict_SetItemString(d, "run_gsl_version", version) != 0){
       fprintf(stderr, "I could not add the run version string to the module dict of pygsl.init!");
       return;
  }

  date = PyString_FromString(DATE " " TIME);
  if(version == NULL){
       fprintf(stderr, "I could not create the date string for pygsl.init!");
       return;
  }
  if(PyDict_SetItemString(d, "compile_date", date) != 0){
       fprintf(stderr, "I could not add the date version string to the module dict of pygsl.init!");
       return;
  }

  if((debuglist = PyList_New(0)) == NULL){
       fprintf(stderr, "Failed to init Debug list!\n");
  }
  /*
   * These functions will be moved to the approbriate modules and the user will
   * have to call them explicitly when needed.
   */

  return;
}

ext: c

max_stars_repo_name: juhnowski/FishingRod

max_issues_repo_name: juhnowski/FishingRod

max_forks_repo_name: juhnowski/FishingRod


text:
/* stable/stable_pdf.c
 *
 * Functions wrappers of GSL routines for random sample generation of
 * alpha-stable random variable.
 *
 * Copyright (C) 2013. Javier Royuela del Val
 *                     Federico Simmross Wattenberg
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; version 3 of the License.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; If not, see <http://www.gnu.org/licenses/>.
 *
 *
 *  Javier Royuela del Val.
 *  E.T.S.I. Telecomunicación
 *  Universidad de Valladolid
 *  Paseo de Belén 15, 47002 Valladolid, Spain.
 *  jroyval@lpi.tel.uva.es
 */
#include "stable.h"
#include <gsl/gsl_randist.h>

void
stable_rnd_seed(StableDist * dist, unsigned long int s)
{
  gsl_rng_set(dist->gslrand, s);
}

inline double
stable_rnd_point(StableDist *dist)
{
  return dist->mu_1 +
         gsl_ran_levy_skew(dist->gslrand, dist->sigma, dist->alpha, dist->beta);
}

void
stable_rnd(StableDist *dist, double *rnd, unsigned int n)
{
  //double *rnd;
  int i;

  //rnd = (double*)malloc(n*sizeof(double));
  if (rnd ==NULL) {
    perror("stable_rnd: NULL output pointer");
    return;
  }

  for(i=0;i<n;i++)
    {
      rnd[i]=stable_rnd_point(dist);
    }
  return;
}

max_stars_repo_name: TantraLabs/gonum

max_issues_repo_name: TantraLabs/gonum

max_forks_repo_name: TantraLabs/gonum


text:
//
// Copyright (C) Microsoft Corporation. All rights reserved.
//

#pragma once

#include <intsafe.h>
#include <gsl/gsl>

namespace arcana
{
    void set_thread_name(DWORD threadId, gsl::czstring<> threadName);
}

ext: h


text:
/* vector/gsl_vector_complex_double.h
 *
 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 Gerard Jungman, Brian Gough
 *
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 3 of the License, or (at
 * your option) any later version.
 *
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
 */

#ifndef __GSL_VECTOR_COMPLEX_DOUBLE_H__
#define __GSL_VECTOR_COMPLEX_DOUBLE_H__

#if !defined( GSL_FUN )
#  if !defined( GSL_DLL )
#    define GSL_FUN extern
#  elif defined( BUILD_GSL_DLL )
#    define GSL_FUN extern __declspec(dllexport)
#  else
#    define GSL_FUN extern __declspec(dllimport)
#  endif
#endif

#include <stdlib.h>
#include <gsl/gsl_types.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_check_range.h>
#include <gsl/gsl_vector_double.h>
#include <gsl/gsl_vector_complex.h>
#include <gsl/gsl_block_complex_double.h>

#undef __BEGIN_DECLS
#undef __END_DECLS
#ifdef __cplusplus
# define __BEGIN_DECLS extern "C" {
# define __END_DECLS }
#else
# define __BEGIN_DECLS /* empty */
# define __END_DECLS /* empty */
#endif

__BEGIN_DECLS

typedef struct
{
  size_t size;
  size_t stride;
  double *data;
  gsl_block_complex *block;
  int owner;
} gsl_vector_complex;

typedef struct
{
  gsl_vector_complex vector;
} _gsl_vector_complex_view;

typedef _gsl_vector_complex_view gsl_vector_complex_view;

typedef struct
{
  gsl_vector_complex vector;
} _gsl_vector_complex_const_view;

typedef const _gsl_vector_complex_const_view gsl_vector_complex_const_view;

/* Allocation */

GSL_FUN gsl_vector_complex *gsl_vector_complex_alloc (const size_t n);
GSL_FUN gsl_vector_complex *gsl_vector_complex_calloc (const size_t n);

GSL_FUN gsl_vector_complex *
gsl_vector_complex_alloc_from_block (gsl_block_complex * b,
                                           const size_t offset,
                                           const size_t n,
                                           const size_t stride);

GSL_FUN gsl_vector_complex *
gsl_vector_complex_alloc_from_vector (gsl_vector_complex * v,
                                             const size_t offset,
                                             const size_t n,
                                             const size_t stride);

GSL_FUN void gsl_vector_complex_free (gsl_vector_complex * v);

/* Views */

GSL_FUN _gsl_vector_complex_view
gsl_vector_complex_view_array (double *base,
                                     size_t n);

GSL_FUN _gsl_vector_complex_view
gsl_vector_complex_view_array_with_stride (double *base,
                                                 size_t stride,
                                                 size_t n);

GSL_FUN _gsl_vector_complex_const_view
gsl_vector_complex_const_view_array (const double *base,
                                           size_t n);

GSL_FUN _gsl_vector_complex_const_view
gsl_vector_complex_const_view_array_with_stride (const double *base,
                                                       size_t stride,
                                                       size_t n);

GSL_FUN _gsl_vector_complex_view
gsl_vector_complex_subvector (gsl_vector_complex *base,
                                         size_t i,
                                         size_t n);

GSL_FUN _gsl_vector_complex_view
gsl_vector_complex_subvector_with_stride (gsl_vector_complex *v,
                                                size_t i,
                                                size_t stride,
                                                size_t n);

GSL_FUN _gsl_vector_complex_const_view
gsl_vector_complex_const_subvector (const gsl_vector_complex *base,
                                               size_t i,
                                               size_t n);

GSL_FUN _gsl_vector_complex_const_view
gsl_vector_complex_const_subvector_with_stride (const gsl_vector_complex *v,
                                                      size_t i,
                                                      size_t stride,
                                                      size_t n);

GSL_FUN _gsl_vector_view
gsl_vector_complex_real (gsl_vector_complex *v);

GSL_FUN _gsl_vector_view
gsl_vector_complex_imag (gsl_vector_complex *v);

GSL_FUN _gsl_vector_const_view
gsl_vector_complex_const_real (const gsl_vector_complex *v);

GSL_FUN _gsl_vector_const_view
gsl_vector_complex_const_imag (const gsl_vector_complex *v);

/* Operations */

GSL_FUN void gsl_vector_complex_set_zero (gsl_vector_complex * v);
GSL_FUN void gsl_vector_complex_set_all (gsl_vector_complex * v,
                                       gsl_complex z);
GSL_FUN int gsl_vector_complex_set_basis (gsl_vector_complex * v, size_t i);

GSL_FUN int gsl_vector_complex_fread (FILE * stream,
                                    gsl_vector_complex * v);
GSL_FUN int gsl_vector_complex_fwrite (FILE * stream,
                                     const gsl_vector_complex * v);
GSL_FUN int gsl_vector_complex_fscanf (FILE * stream,
                                     gsl_vector_complex * v);
GSL_FUN int gsl_vector_complex_fprintf (FILE * stream,
                                      const gsl_vector_complex * v,
                                      const char *format);

GSL_FUN int gsl_vector_complex_memcpy (gsl_vector_complex * dest, const gsl_vector_complex * src);

GSL_FUN int gsl_vector_complex_reverse (gsl_vector_complex * v);

GSL_FUN int gsl_vector_complex_swap (gsl_vector_complex * v, gsl_vector_complex * w);
GSL_FUN int gsl_vector_complex_swap_elements (gsl_vector_complex * v, const size_t i, const size_t j);

GSL_FUN int gsl_vector_complex_equal (const gsl_vector_complex * u,
                                    const gsl_vector_complex * v);

GSL_FUN int gsl_vector_complex_isnull (const gsl_vector_complex * v);
GSL_FUN int gsl_vector_complex_ispos (const gsl_vector_complex * v);
GSL_FUN int gsl_vector_complex_isneg (const gsl_vector_complex * v);
GSL_FUN int gsl_vector_complex_isnonneg (const gsl_vector_complex * v);

GSL_FUN int gsl_vector_complex_add (gsl_vector_complex * a, const gsl_vector_complex * b);
GSL_FUN int gsl_vector_complex_sub (gsl_vector_complex * a, const gsl_vector_complex * b);
GSL_FUN int gsl_vector_complex_mul (gsl_vector_complex * a, const gsl_vector_complex * b);
GSL_FUN int gsl_vector_complex_div (gsl_vector_complex * a, const gsl_vector_complex * b);
GSL_FUN int gsl_vector_complex_scale (gsl_vector_complex * a, const gsl_complex x);
GSL_FUN int gsl_vector_complex_add_constant (gsl_vector_complex * a, const gsl_complex x);

GSL_FUN INLINE_DECL gsl_complex gsl_vector_complex_get (const gsl_vector_complex * v, const size_t i);
GSL_FUN INLINE_DECL void gsl_vector_complex_set (gsl_vector_complex * v, const size_t i, gsl_complex z);
GSL_FUN INLINE_DECL gsl_complex *gsl_vector_complex_ptr (gsl_vector_complex * v, const size_t i);
GSL_FUN INLINE_DECL const gsl_complex *gsl_vector_complex_const_ptr (const gsl_vector_complex * v, const size_t i);

#ifdef HAVE_INLINE

INLINE_FUN
gsl_complex
gsl_vector_complex_get (const gsl_vector_complex * v,
                              const size_t i)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(i >= v->size))
    {
      gsl_complex zero = {{0, 0}};
      GSL_ERROR_VAL ("index out of range", GSL_EINVAL, zero);
    }
#endif
  return *GSL_COMPLEX_AT (v, i);
}

INLINE_FUN
void
gsl_vector_complex_set (gsl_vector_complex * v,
                              const size_t i, gsl_complex z)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(i >= v->size))
    {
      GSL_ERROR_VOID ("index out of range", GSL_EINVAL);
    }
#endif
  *GSL_COMPLEX_AT (v, i) = z;
}

INLINE_FUN
gsl_complex *
gsl_vector_complex_ptr (gsl_vector_complex * v,
                              const size_t i)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(i >= v->size))
    {
      GSL_ERROR_NULL ("index out of range", GSL_EINVAL);
    }
#endif
  return GSL_COMPLEX_AT (v, i);
}

INLINE_FUN
const gsl_complex *
gsl_vector_complex_const_ptr (const gsl_vector_complex * v,
                                    const size_t i)
{
#if GSL_RANGE_CHECK
  if (GSL_RANGE_COND(i >= v->size))
    {
      GSL_ERROR_NULL ("index out of range", GSL_EINVAL);
    }
#endif
  return GSL_COMPLEX_AT (v, i);
}

#endif /* HAVE_INLINE */

__END_DECLS

#endif /* __GSL_VECTOR_COMPLEX_DOUBLE_H__ */


text:
/*
 * rw.c -- Loesung des Randwertproblems mit der GSL
 *
 * (c) 2016 Prof Dr Andreas Mueller, Hochschule Rapperswil
 */
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_odeiv2.h>

static double	g = 9.81;

static long	fcounter = 0;

static int	f(double x, const double y[], double f[], void *params) {
	fcounter++;
	double	g = *(double *)params;
	f[0] =  y[2];
	f[1] =  y[3];
	f[2] =     0;
	f[3] =    -g;
	return GSL_SUCCESS;
}

static long	dfcounter = 0;

static int	df(double x, const double y[], double *dfdy, double dfdx[],
			void *params) {
	dfcounter++;

	dfdy[0 * 4 + 0] = 0;
	dfdy[1 * 4 + 0] = 0;
	dfdy[2 * 4 + 0] = 0;
	dfdy[3 * 4 + 0] = 0;

	dfdy[0 * 4 + 1] = 0;
	dfdy[1 * 4 + 1] = 0;
	dfdy[2 * 4 + 1] = 0;
	dfdy[3 * 4 + 1] = 0;

	dfdy[0 * 4 + 2] = 1;
	dfdy[1 * 4 + 2] = 0;
	dfdy[2 * 4 + 2] = 0;
	dfdy[3 * 4 + 2] = 0;

	dfdy[0 * 4 + 3] = 0;
	dfdy[1 * 4 + 4] = 1;
	dfdy[2 * 4 + 4] = 0;
	dfdy[3 * 4 + 4] = 0;

	dfdx[0] = 0;
	dfdx[1] = 0;
	dfdx[2] = 0;
	dfdx[3] = 0;

	return GSL_SUCCESS;
}

static long	Fcounter = 0;

static int	F(double x, const double J[], double F[], void *params) {
	Fcounter++;
	for (int i = 0; i < 2; i++) {
		for (int j = 0; j < 4; j++) {
			F[i * 4 + j] = J[(i + 2) * 4 + j];
		}
	}
	for (int i = 2; i < 4; i++) {
		for (int j = 0; j < 4; j++) {
			F[i * 4 + j] = 0;
		}
	}
	return GSL_SUCCESS;
}

int	main(int argc, char *argv[]) {
	gsl_odeiv2_system	system = { f, df, 4, &g };
	gsl_odeiv2_system	Jsystem = { F, NULL, 16, NULL };

	gsl_odeiv2_driver	*driver
		= gsl_odeiv2_driver_alloc_y_new(&system, gsl_odeiv2_step_rk8pd,
			1e-6, 1e-6, 0.0);
	gsl_odeiv2_driver	*Jdriver
		= gsl_odeiv2_driver_alloc_y_new(&Jsystem, gsl_odeiv2_step_rk8pd,
			1e-6, 1e-6, 0.0);

	double	t = 0.0;
	double	vx = 8;
	double	vy = 7;
	double	tnext = 20 / vx;
	double	delta;
	int	counter = 0;
	printf(" n     v_y        t        x        y        dy/dv_y        vynew          delta\n");
	do {

		// compute the solution up to x = 1.5
		t = 0;
		double	y[4] = { 0.0, 0.0, vx, vy };
		int	status = gsl_odeiv2_driver_apply(driver, &t, tnext, y);
		if (status != GSL_SUCCESS) {
			fprintf(stderr, "error: return value = %d\n", status);
		}

		// compute the jacobi matrix up to x = 1.5
		t = 0;
		double	J[16] = {	1.0, 0.0, 0.0, 0.0,
					0.0, 1.0, 0.0, 0.0,
					0.0, 0.0, 1.0, 0.0,
					0.0, 0.0, 0.0, 1.0 };
		status = gsl_odeiv2_driver_apply(Jdriver, &t, tnext, J);
		if (status != GSL_SUCCESS) {
			fprintf(stderr, "error: return value = %d\n", status);
		}

		// the derivative is J_{24}
		double	 derivative = J[1 * 4 + 3];

		// compute the newton correction
		double	vynew = vy - y[1] / derivative;
		delta = vy - vynew;
		printf("%2d%8.4f %8.4f %8.4f %10.6f %12.8f %12.8f %14.10f\n",
			 counter, vy, t, y[0], y[1], derivative, vynew, delta);
		vy = vynew;
		counter++;
	} while ((fabs(delta) > 1e-6) && (counter <= 50));
	printf("calls to f: %ld, to df: %ld, to F: %ld\n",
		fcounter, dfcounter, Fcounter);
	return EXIT_SUCCESS;
}

max_stars_repo_name: MatthiasRubin/SeminarDGL

max_issues_repo_name: MatthiasRubin/SeminarDGL

max_forks_repo_name: MatthiasRubin/SeminarDGL


text:
#ifndef GSLUtility_h
#define GSLUtility_h
/** This is free and unencumbered software released into the public domain.
The authors of ISIS do not claim copyright on the contents of this file.
For more details about the LICENSE terms and the AUTHORS, you will
find files of those names at the top level of this repository. **/

/* SPDX-License-Identifier: CC0-1.0 */

#include <iostream>
#include <sstream>
#include <string>
#include <vector>
#include <algorithm>

#include <QString>

#include <tnt/tnt_array1d.h>
#include <tnt/tnt_array1d_utils.h>
#include <tnt/tnt_array2d.h>
#include <tnt/tnt_array2d_utils.h>

#include <gsl/gsl_vector.h>
#include <gsl/gsl_matrix.h>

//  Some GSL optimization on by default, off if DEBUG or SAFE_GSL is defined
#ifndef DEBUG
#ifndef SAFE_GSL
#define GSL_RANGE_CHECK_OFF 1
#endif
#endif

#include "IException.h"

namespace Isis {
  namespace GSL {

    /**
     * @brief GSLUtility Provides top level interface to the GNU GSL
     *
     * Provides GSL setup and interface utilities.  This object is provided for
     * convenience of GSL vector and matrix manipulation as well as better
     * management of GSL error handling.
     *
     * Without setting up GSL error handling, the GSL will abort when certain
     * errors occur.  This singleton object, an object where there is never more
     * than one instance, an error handler is established that captures GSL errors
     * and formats them into ISIS exceptions.
     *
     * There are many convenience methods provided for manipulation of GSL vectors
     * and matrixs.  Motivation for this is to address element access and
     * efficient parameter and copy mechanisms (provided by the TNT library).
     *
     * There are some compile options on by default that help optimize the GSL.
     * When the compile time DEBUG macro is set, range checking is turned on.
     * Inline functions are also turned on by default unless the DEBUG macro is
     * set. In addition, an additional compile time macro called SAFE_GSL is
     * provided to emulate the DEBUG behavior, but does not invoke additional
     * DEBUG behavior/side effects.
     *
     * See http://www.gnu.org/software/gsl/ for additional details on the GNU
     * Scientific Library.
     *
     * @ingroup Utility
     * @author 2008-05-06 Kris Becker
     * @internal
     *   @history 2009-08-20 Kris Becker Completed documentation
     */
    class GSLUtility {
      public:
        typedef TNT::Array1D<double> GSLVector;
        typedef TNT::Array2D<double> GSLMatrix;

        static GSLUtility *getInstance();

        /** Tests if status is success */
        inline bool success(int status) const {
          return (status == GSL_SUCCESS);
        }

        /**
         * @brief Returns GSL specific error text
         *
         * @param gsl_errno GSL error number
         * @return QString Textual context of GSL error
         */
        inline QString status(int gsl_errno) const {
          return (QString(gsl_strerror(gsl_errno)));
        }

        void check(int gsl_status, const char *src = __FILE__, int line = __LINE__)
        const;

        size_t Rows(const gsl_matrix *m) const;
        size_t Columns(const gsl_matrix *m) const;

        size_t Rows(const GSLMatrix &m) const;
        size_t Columns(const GSLMatrix &m) const;

        size_t size(const gsl_vector *v) const;
        size_t size(const gsl_matrix *m) const;

        gsl_vector *vector(size_t n, bool zero = false) const;
        gsl_matrix *matrix(size_t n1, size_t n2, bool zero = false) const;
        gsl_matrix *identity(size_t n1, size_t n2) const;
        void setIdentity(gsl_matrix *m) const;

        void free(gsl_vector *v) const;
        void free(gsl_matrix *m) const;

        GSLVector gslToGSL(const gsl_vector *v) const;
        GSLMatrix gslToGSL(const gsl_matrix *m) const;
        gsl_vector *GSLTogsl(const GSLVector &v, gsl_vector *gv = 0) const;
        gsl_matrix *GSLTogsl(const GSLMatrix &m, gsl_matrix *gm = 0) const;

      private:
        //  Private Constructor/Destructor makes this a singleton
        GSLUtility();
        ~GSLUtility() { }

        static GSLUtility *_instance;  //!< Singleton self-reference pointer

        static void handler(const char *reason, const char *file, int line,
                            int gsl_errno);

    };

  } // namespace GSL
}     // namespace Isis
#endif


text:
#include <stdio.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_matrix.h>
#include <time.h>

#define N_OF_RUNS 10

#define RUSAGE_DIFF(r1, r2) \
    ((double) ((r1).tv_sec - (r2).tv_sec) + (double) ((r1).tv_nsec - (r2).tv_nsec)/(double)1E9)

gsl_matrix* random_matrix(size_t);
inline gsl_matrix* zero_matrix(const size_t size);

double* naive_test(size_t);
double* better_test(size_t);
double* dgemm_test(size_t);

void test_suite(size_t, size_t, size_t, const char*);

int main() {
    srand48((unsigned int)time(NULL));
    const char filename[] = "results_unoptimized.csv";
    test_suite(100, 1000, 100, filename);
    return 0;
}

void test_suite(const size_t from, const size_t to, const size_t step, const char* filename){
    FILE* file = fopen(filename, "w");
    if(file == NULL){
        printf("Error opening file, exiting \r\n");
        exit(EXIT_FAILURE);
    }
    fprintf(file, "METHOD,SIZE,TIME\r\n");

    double* naive_result, * dgemm_result, * better_result;

    for(size_t i = from; i <= to; i += step){
        for(int j = 0; j < N_OF_RUNS; j++) {
            naive_result = naive_test(i);
            better_result = better_test(i);
            dgemm_result = dgemm_test(i);
            fprintf(file, "Naive,%d,%lf\r\n", (int) i, *naive_result);
            fprintf(file, "Better,%d,%lf\r\n", (int) i, *better_result);
            fprintf(file, "BLAS,%d,%lf\r\n", (int) i, *dgemm_result);
            free(naive_result);
            free(better_result);
            free(dgemm_result);
        }
    }

    fclose(file);
}

gsl_matrix* random_matrix(const size_t size) {
    gsl_matrix* a = gsl_matrix_calloc(size, size);
    for(int i = 0; i < size; i++)
        for(int j = 0; j < size; j++)
            gsl_matrix_set(a, (size_t)i, (size_t)j, drand48()*10);
    return a;
}

gsl_matrix* zero_matrix(const size_t size) {
    gsl_matrix* a = gsl_matrix_calloc(size, size);
    return a;
}

double** random_double_matrix(const size_t size) {
    double** a = calloc(size, sizeof(double*));
    for(int i = 0; i < size; i++)
        a[i] = calloc(size, sizeof(double));

    for(int i = 0; i < size; i++)
        for(int j = 0; j < size; j++)
            a[i][j] = drand48() * 10;
    return a;
}

double** zero_double_matrix(const size_t size) {
    double** a = calloc(size, sizeof(double*));
    for(int i = 0; i < size; i++)
        a[i] = calloc(size, sizeof(double));
    return a;
}

void free_double_matrix(double** matrix, const size_t size){
    for(int i = 0; i < size; i++)
        free(matrix[i]);
    free(matrix);
}

double* naive_test(const size_t size) {
    double** m_1 = random_double_matrix(size);
    double** m_2 = random_double_matrix(size);
    double** m_3 = zero_double_matrix(size);

    printf("Multiplying matrices of size %d using naive method\r\n", (int)size);

    struct timespec real_start, real_end;
    clock_gettime(CLOCK_REALTIME, &real_start);

    for(int k = 0; k < size; k++)
        for(int j = 0; j < size; j++)
            for(int i = 0; i < size; i++)
                m_3[i][j] += m_1[i][k] * m_2[k][j];

    clock_gettime(CLOCK_REALTIME, &real_end);

    double* time = malloc(sizeof(double));
    *time = RUSAGE_DIFF(real_end, real_start);
    printf("Real time = %lf \r\n", *time);
    free_double_matrix(m_1, size); free_double_matrix(m_2, size); free_double_matrix(m_3, size);
    return time;
}

double* better_test(const size_t size) {
    double** m_1 = random_double_matrix(size);
    double** m_2 = random_double_matrix(size);
    double** m_3 = zero_double_matrix(size);

    printf("Multiplying matrices of size %d using better method\r\n", (int)size);

    struct timespec real_start, real_end;
    clock_gettime(CLOCK_REALTIME, &real_start);

    for(int i = 0; i < size; i++)
        for(int j = 0; j < size; j++)
            for(int k = 0; k < size; k++)
                m_3[i][j] += m_1[i][k] * m_2[k][j];

    clock_gettime(CLOCK_REALTIME, &real_end);

    double* time = malloc(sizeof(double));
    *time = RUSAGE_DIFF(real_end, real_start);
    printf("Real time = %lf \r\n", *time);
    free_double_matrix(m_1, size); free_double_matrix(m_2, size); free_double_matrix(m_3, size);
    return time;
}

double* dgemm_test(const size_t size) {
    gsl_matrix* m_1 = random_matrix(size);
    gsl_matrix* m_2 = random_matrix(size);
    gsl_matrix* m_3 = zero_matrix(size);

    printf("Multiplying matrices of size %d using BLAS\r\n", (int)size);

    struct timespec real_start, real_end;
    clock_gettime(CLOCK_REALTIME, &real_start);

    gsl_blas_dgemm (CblasNoTrans, CblasNoTrans, 1.0, m_1, m_2, 0.0, m_3);

    clock_gettime(CLOCK_REALTIME, &real_end);

    double* time = malloc(sizeof(double));
    *time = RUSAGE_DIFF(real_end, real_start);
    printf("Real time = %lf \r\n", *time);
    gsl_matrix_free(m_1); gsl_matrix_free(m_2); gsl_matrix_free(m_3);
    return time;
}


text:
/*
* -----------------------------------------------------------------
*  ell_lib.h
*  Ellipsoid Library
*  Version: 2.0
*  Last Update: July 22, 2020
*
*  Programmer: Americo Barbosa da Cunha Junior
*              americo.cunhajr@gmail.com
* -----------------------------------------------------------------
*  Copyright (c) 2010-2019, Americo Barbosa da Cunha Junior
*  All rights reserved.
* -----------------------------------------------------------------
*  This is the header file for BST_LIB module, a computational
*  library to work with n-dimensional ellipsoids.
* -----------------------------------------------------------------
*/

#ifndef __ELL_LIB_H__
#define __ELL_LIB_H__

#include <stdlib.h>
#include <gsl/gsl_linalg.h>

#undef  ELL_TOL
#define ELL_TOL 1.0e-7

#undef  ELL_TRUE
#define ELL_TRUE 1

#undef  ELL_FALSE
#define ELL_FALSE 0

/*
*------------------------------------------------------------
* function prototypes
*------------------------------------------------------------
*/

void ell_psd2eig(gsl_matrix *B,
                 gsl_matrix *V,
                 gsl_vector *lam);

void ell_psd2chol(gsl_matrix *B,
                  gsl_matrix *L);

void ell_eig2chol(gsl_matrix *V,
                  gsl_vector *sig,
                  gsl_matrix *L);

unsigned int ell_pt_in(gsl_vector *x,
                       gsl_vector *c,
                       gsl_matrix *L);

void ell_map2uhs(gsl_vector *x,
                 gsl_vector *c,
                 gsl_matrix *L,
                 gsl_vector *y);

void ell_pt_modify(gsl_vector *p,
                   gsl_vector *c,
                   gsl_matrix *L);

#endif /* __ELL_LIB_H__ */

max_stars_repo_name: americocunhajr/CRFlowLib

max_issues_repo_name: americocunhajr/CRFlowLib

max_forks_repo_name: americocunhajr/CRFlowLib


text:
#include <gsl/gsl_integration.h>
#include <pygsl/error_helpers.h>
#include <pygsl/function_helpers.h>
#include <setjmp.h>

PyObject *module = NULL;
typedef struct{
	PyObject * callback;
	PyObject * args;
	jmp_buf  buffer;
}pygsl_function_args;

static double
pygsl_function_callback(double x, void *p)
{
	double value;
	int flag;
	pygsl_diff_args *pargs = NULL;

	pargs = (pygsl_diff_args *) p;

	assert(pargs->callback);
	assert(pargs->args);
	flag = PyGSL_function_wrap_helper(x, &value, NULL, pargs->callback,
					  pargs->args, (char *)__FUNCTION__);

	if(GSL_SUCCESS != flag){
		longjmp(pargs->buffer, flag);
		return gsl_nan();
	}
	return value;
}

PyObject *
integrate(PyObject *self, PyObject *pyargs)
{
     PyObject *callback = NULL, *args = NULL;
     pygsl_function_args cb_args;
     gsl_function f;
     int flag;
     double a, b, epsabs, epsrel, result, abserr;
     size_t neval;

     if(!PyArg_ParseTuple(args, "Odddd|O", &callback, &a, &b, &epsabs, &epsrel, &args))
	  goto fail;

     if(!Py_Callable(callback)){
	  pygsl_error("First argument must be a callable object!", GSL_EINVAL);
	  goto fail;
     }

     if(args == NULL){
	  args = Py_NONE;
     }

     f.f = pygsl_function_callback;
     f.params = &cb_args;
     cb_args.args = args;
     cb_args.callback = callback;
     flag = gsl_integration_qng (const *F, a, b, epsabs, epsrel, &result, &abserr, &neval)

     if(PyGSL_ERROR_FLAG(flag) != GSL_SUCCESS)
	  goto fail;

     return Py_BuildValue("(ddi)", result, abserr, neval);

 fail:
     Py_XDECREF(callback);
     Py_XDECREF(args);
     return NULL;
}
void
initintegrate(void)
{
     PyObject *m;

     m=Py_InitModule("integrate", mMethods);
     module = m;
     assert(m);
}