/* part of the shansyn spherical harmonics package, see COPYRIGHT for license */
/* $Id: write_coefficients.c,v 1.7 2004/06/22 00:31:28 becker Exp becker $ */
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include "function_macros.h"
#include "trig_constants.h"
#include "precision.h"
#include "spherical_harmonics_functions.h"
#include "legendre_macros.h"


void write_nonzero_coefficients(COMP_PRECISION **a,
				COMP_PRECISION **b,int lmax,
				COMP_PRECISION sfac,int print_lmax,
				FILE *out, int nexp)
{
  int l,m,i,os;
  if(print_lmax)
    fprintf(out,"%i\n",lmax);
  for(l=0;l <= lmax;l++)
    for(m=0;m <= l;m++){
      os = POSLM(l,m);
      if(m == 0){
	/* don't write the zero B term for m == 0 */
	for(i=0;i < nexp;i++)
	  fprintf(out,ASCII_DATA_FORMAT,*(a[i]+os));
      }else{
	for(i=0;i < nexp;i++)
	  fprintf(out,COEFF_DATA_FORMAT,*(a[i]+os),*(b[i]+os));
      }
      fprintf(out,"\n");
    }

}
void write_coefficients(COMP_PRECISION **a,
			COMP_PRECISION **b,int lmax,
			COMP_PRECISION sfac,int print_lmax,
			FILE *out, int nexp)
{
  int l,m,i,os;
  if(print_lmax)
    fprintf(out,"%i\n",lmax);
  if(sfac != 1.0){
    for(l=0;l<=lmax;l++)
      for(m=0;m<=l;m++){
	os = POSLM(l,m);
	for(i=0;i<nexp;i++){
	  if((*(a[i]+os) != 0)||(*(b[i]+os)!=0))
	    fprintf(out,COEFF_DATA_FORMAT,*(a[i]+os)*sfac,*(b[i]+os)*sfac);
	  else
	    fprintf(out,"              0               0");
	}
	fprintf(out,"\n");
      }
  }else{
    for(l=0;l<=lmax;l++)
      for(m=0;m<=l;m++){
	os = POSLM(l,m);
	for(i=0;i<nexp;i++){
	  if((*(a[i]+os) != 0)||(*(b[i]+os)!=0))
	    fprintf(out,COEFF_DATA_FORMAT,*(a[i]+os),*(b[i]+os));
	  else
	    fprintf(out,"              0               0");
	}
	fprintf(out,"\n");
      }
  }
}
void write_model(struct mod *model,COMP_PRECISION fac,
		 int nexp)
{
  int i,j;
  COMP_PRECISION **a,**b;
  a=(COMP_PRECISION **)malloc(sizeof(COMP_PRECISION *)*nexp);
  b=(COMP_PRECISION **)malloc(sizeof(COMP_PRECISION *)*nexp);

  switch(model[0].type){
  case DISCRETE:{
    printf("%i\n",model[0].n);// number of layers
    for(i=0;i < model[0].n;i++){
      for(j=0;j<nexp;j++){
	a[j] = model[j].a[i];
	b[j] = model[j].b[i];
      }
      printf("%.8e\n",model[0].d[i]);// depth
      // coefficients
      write_coefficients(a,b,model[0].lmax,fac,TRUE,stdout,nexp);
    }
    break;
  }
  case CHEBYSHEV:{
    printf("%i\n",-model[0].n);// use negative order as layer number
    printf("%.8e %.8e\n",model[0].dmin,model[0].dmax);// boundaries
    for(i=0;i<model[0].n;i++){// coefficients
      for(j=0;j<nexp;j++){
	a[j] = model[j].a[i];
	b[j] = model[j].b[i];
      }
      write_coefficients(a,b,model[0].lmax,fac,TRUE,stdout,nexp);
    }
    break;
  }
  case SPLINES:{
    printf("%i\n",model[0].n+1000);// use order + 1000 as layer number
    printf("%.8e %.8e\n",model[0].dmin,model[0].dmax);// boundaries
    for(i=0;i<model[0].n;i++){// coefficients
      for(j=0;j<nexp;j++){
	a[j] = model[j].a[i];
	b[j] = model[j].b[i];
      }
      write_coefficients(a,b,model[0].lmax,fac,TRUE,stdout,nexp);
    }
    break;
  }
  case HC_TYPE:{
    for(i=0;i < model[0].n;i++){
      fprintf(stdout,"%i %i %.8e %i %i %i\n",
	      model[0].lmax,i,model[0].d[i],model[0].n,1,0);
      for(j=0;j<nexp;j++){
	a[j] = model[j].a[i];
	b[j] = model[j].b[i];
      }
      write_coefficients(a,b,model[0].lmax,fac,FALSE,stdout,nexp);
    }
    break;
  }
  default:{
    fprintf(stderr,"write_model: type %i undefined\n",
	    model[0].type);
    exit(-1);
  }}
  free(a);free(b);
}
