/*
 This C routine outputs model solution data cube with the correct format for the SHINE validation study

 The arguments to this routine are

  nr                integer                                   number of radial grid points
  nth               integer                                   number of latitudinal grid points
  nph               integer                                   number of longitudinal grid points
  r                 float vector    dimension nr              radial grid point coords 
  th                float vector    dimension nth             latitudinal grid point coords 
  ph                float vector    dimension nph             longitudinal grid point coords 
  br                float vector    dimension nr*nth*nph      vector storing radial component of magnetic field
  bth               float vector    dimension nr*nth*nph      vector storing latitudinal component of magnetic field
  bph                float vector   dimension nr*nth*nph      vector storing longitudinal component of magnetic field
  dens              float vector    dimension nr*nth*nph      vector storing density at grid points
  temp              float vector    dimension nr*nth*nph      vector storing temperature at grid points
  pres              float vector    dimension nr*nth*nph      vector storing pressure at grid points
  units_descriptor  character string                          string describing units of variables
  filename          character string                          string containing the output filename


  The 3D data in br,bth,bph,dens,temp and pres is passed into this routine in a flar 1D vector with
  the radial coordinate index changing fastest, followed by the latitudinal index and finally the longitudinal
  index changing slowest.


  Written : Peter MacNeice                March 2012
*/

#include <stdio.h>
#include <string.h>
#include <stdlib.h>


void shine_valid_format_c( 
     int nr, int nth, int nph,
     float *r, float *th, float *ph,
     float *br, float *bth, float *bph, 
     float *dens, float *temp, float *pres,
     char *units_descriptor,
     char *filename
                         )

{

    FILE * fp;

    int i, j, k, l;

    fp = fopen(filename,"w");


    fprintf(fp, "%s\n",units_descriptor);
    fprintf(fp, "%d %d %d\n",nr,nth,nph);

    for ( i=0; i<nr; i++) {
      fprintf(fp, "%e\n",r[i]);
    }

    for ( i=0; i<nth; i++) {
      fprintf(fp, "%e\n",th[i]);
    }

    for ( i=0; i<nph; i++) {
      fprintf(fp, "%e\n",ph[i]);
    }


    for ( k=0; k<nph; k++) {
      for ( j=0; j<nth; j++) {
        for ( i=0; i<nr;  i++)  {
          l = i + nr*(j + nth*k);
          fprintf(fp, "%e %e %e\n",br[l], bth[l], %bph[l]);
        }
      }
    }

    for ( k=0; k<nph; k++) {
      for ( j=0; j<nth; j++) {
        for ( i=0; i<nr; i++)  {
          l = i + nr*(j + nth*k);
          fprintf(fp, "%e %e %e\n",dens[l], %temp[l], %pres[l]);
        }
      }
    }

    fclose(fp);

}