#!/usr/bin/env /home/paco/hacking/cint/runme 
-lgsl -lgslcblas -lm

#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <stdarg.h>
#include <gsl/gsl_sf_bessel.h>

// The multidimensional array struct which everything is.
typedef struct structarray {
	int d;		// number of dimensions
	int *s;		// shape vector	
	double *b;	// block of data
} sa;

// Fill an array with a constant value.
sa fc(sa p, double val)  
{
    int dim;
    dim=p.d;
    int lcv;
    int totnum=1;
    for (lcv=0;lcv<dim;lcv++){
	totnum*=p.s[lcv];
    }    
    int lcv2;
    for (lcv2=0;lcv2<totnum;lcv2++)
    	p.b[lcv2]= val;
    return p;
}

// Fill an array with an incremental range.
sa fr(sa p)  
{
    int dim;
    dim=p.d;
    int lcv;
    int totnum=1;
    for (lcv=0;lcv<dim;lcv++){
	totnum*=p.s[lcv];
    }    
    int lcv2;
    for (lcv2=0;lcv2<totnum;lcv2++)
    	p.b[lcv2]= (double) lcv2;
    return p;
}

//Create a new array given the shape 
sa cn(int dim, ...)  
{
    sa p;
    p.d =dim;
    p.s=  (int *) malloc(dim*sizeof(int)); 
    va_list ap;
    va_start(ap,dim);
    int lcv;
    int totnum=1;
    for (lcv=0;lcv<dim;lcv++){
    	p.s[lcv]= va_arg(ap,int);
	totnum*=p.s[lcv];
    }    
    va_end(ap);
    p.b=  (double *) malloc(totnum*sizeof(double)); 
    return p;
}

//Create a new array with the same shape as an old one.
sa cs(sa a)  
{
    sa p;
    p.d =a.d;
    p.s=  (int *) malloc(p.d*sizeof(int)); 
    int lcv;
    int totnum=1;
    for (lcv=0;lcv<p.d;lcv++){
    	p.s[lcv]=a.s[lcv];
	totnum*=a.s[lcv];
    }    
    p.b=  (double *) malloc(totnum*sizeof(double)); 
    return p;
}

// Helper function to pretty-print.
void pb(int *elptr,int lev,int *p, sa a){
	if (lev>0){
		printf("[");
	for (p[lev-1]=0;p[lev-1]<a.s[lev-1];p[lev-1]++)
		pb(elptr,lev-1,p,a);
		printf("]\n");
	}
	else {
		printf("%5.2f ",a.b[*elptr]);
		(*elptr)++;
	}
}

// Pretty-print an array.
void pa (sa a){
	printf("Dimensions: %d\n",a.d);
	printf("Shape: ");
	int lcv;
	for (lcv=0;lcv<a.d;lcv++){
	    printf("%d ",a.s[lcv]);
	}
	printf("\n");
	int *p= (int *)malloc(a.d*sizeof(int));
	int elptr=0;
	pb(&elptr,a.d,p,a);
}

// Print the shape of an array.
void ps (sa a){
	printf("Dimensions: %d\n",a.d);
	printf("Shape: ");
	int lcv;
	for (lcv=0;lcv<a.d;lcv++){
	    printf("%d ",a.s[lcv]);
	}
	printf("\n");
}

// Create an array by Reading a file
sa cr(char *s);

// Apply a unary operator function to a whole array. 
sa du(sa a,double (*uop)(double)){
	sa p;
	p=cs(a);
	int lcv,lcv2,totnum=1;
	for (lcv=0;lcv<p.d;lcv++){
		totnum*=p.s[lcv];
	}
	for (lcv2=0;lcv2<totnum;lcv2++)
		p.b[lcv2]=(*uop)(a.b[lcv2]);
	return p;
}

//Main function to test things.
int main(int argc, char* argv[]) {
  // Allocate a vector of pointers to arrays.
  sa a[100];
  // Create an array with one dimension, 10 elements, fill it with a range.
  a[1]=fr(cn(1,10));
  // Print the shape.
  ps(a[1]);
  // Print it.
  pa(a[1]);
  // Apply the unary zero-order bessel function to an array.
  a[2]=du(a[1],gsl_sf_bessel_J0);
  // Print it.
  pa(a[2]);
  // Create an array of a constant value.
  a[3]=fc(cn(1,10),6.5);
  // Print it.
  pa(a[3]);
  // Fill with a bunch of values, fv(sa, 1,2,3,4);
  // Fill with a string, fs(sa, "1 2 3 4");
  // Reshape an array, preserving total number of elements.
  // Slice an array, preserving number of dimensions.
  // Concatenate two arrays, checking that dimensions match.
  // Apply binary operators to two arrays. (check matching size first)
  // Plot sin(x) data using gnuplot. 
  // load a winspec data file, the 0 vector will be lambda, 1.. will be data.
  // Fit using gsl fitting routines.
  // Matrix-Multiply two arrays together in a way that is compatible with gsl
  // routines (either use them or match the result).
  // Find the index of maximum value along an axis. max(sa, axis).
  // Find the index of minimum value along an axis. min(sa, axis).
  //  (for the above two, return an array missing the axis requested)
  printf("Hello, World\n");
  // Return 0 for no error. 
  return 0;
}