/* compiles by `gcc -lm dust.c -o dust' */

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

#define MYVERSION 1.04
#define Pi 3.141592653589793
#define NUM_OF_DEF 6
#define MASS_COMPONENT 5
#define K_MAX 10000000
#define MAX_NUM_OF_COMP 10
#define MAX_PARSE 10
#define MAX_LENGTH 256
#define LAMBDA_DEFAULT 550.
#define MAX_LENGTH_FLINE 80

#define ICE_A  7.87e39
#define ICE_al 0.48175
#define ICE_B  8100.
#define ICE_bl 6700.

#define float double
#define printff ZERO_MESSAGES ? : printf

float SIGMA_RANGE_MINUS=3.;
float SIGMA_RANGE_PLUS=7.;
int MAX_STEPS=600;
unsigned char INCLUDE_ICE=0;
unsigned char NO_MESSAGES=0;
unsigned char ZERO_MESSAGES=0;
unsigned char ALLOW=0;
unsigned char USE_INDEX=0;
unsigned char VERBAL=0;
unsigned char DEBUG=0;
char fast_flag=0;
float ang=0;

int line_input(char *str)
{
  int i;
  char ch;
  if(!NO_MESSAGES) printff(">  ");
  for(i=0;i<MAX_LENGTH-1;i++){
    ch=getchar();
    if(ch!='\n') str[i]=ch;
    else break;
  }
  if(i==MAX_LENGTH-1){
    printff("\
Warning: your line exceeded %i characters; it will be cut to fit the buffer\n", MAX_LENGTH-1);
    while(getchar()!='\n');
  }
  str[i]='\0';
  return i;
}


int input(float components[MAX_NUM_OF_COMP][NUM_OF_DEF], int num)
{
  char *help="This is a parser for the scattering coefficients program, v%g\n\
By default, you need to supply %i numbers corresponding to mass densities of\n\
each of the %i components in each input line. If you make an error, you will\n\
be warned that your input is rejected, and given the opportunity to correct\n\
the error. You can always quit the parser by entering keywords `exit' or\n\
`quit'. This results in an exit code  being passed along to the calling\n\
function, which will most likely exit the program. Use `#' for comments.\n\
All keywords available are:\n\
\n\
    `help'  print this help message\n\
    `exit' = `quit' quit the parser\n\
    `show'  show current settings\n\
    `set'   change current parsing settings\n\
    `susp'  turn suspended mode on/off\n\
    `exec'  execute now\n\
    `allow'     | allow/disallow execution\n\
    `disallow'  |  of the system commands\n\
\n";
  char *onoff[2]={"off", "on"};
  char *allow_disallow[2]={"disallow", "allow"};

  int i, j;
  unsigned char flag;
  static unsigned char state=0, suspend=0;
  static int lnum;
  int lnum_aux;
  static int mask[MAX_PARSE];
  int mask_aux[MAX_PARSE];
  static float mass[MAX_NUM_OF_COMP];
  float mass_aux[MAX_NUM_OF_COMP];
  static int mass_comp=MASS_COMPONENT;
  static int num_mass=0;
  int length;
  int iaux;
  int position;
  char str[MAX_LENGTH], key[MAX_LENGTH];
  char *saux;

  if(!VERBAL) return 0;
  if(state==0){
    state=1;
    printff("You need to provide %i mass densities in each line\n\
exit quits, type help for help\n", num);
    lnum=num;
    for(i=0;i<num;i++) mask[i]=i+1;
    for(i=num;i<MAX_PARSE;i++) mask[i]=0;
    for(i=0;i<num;i++) mass[i]=0;
  }
  else if(state==1) state=2;

  while(1){
    length=line_input(str);
    if(sscanf(str, "%s%n", key, &position)!=1) continue;

    if(!strcmp(key, "help")) printff(help, MYVERSION, num, num);
    else if(!(strcmp(key, "exit")&&strcmp(key, "quit"))) return 1;

    else if(!strcmp(key, "allow")){
      if(sscanf(str+position, "%s%n", key, &position)!=1){
	ALLOW=1;
	printff("Execution of external (system) commands will be allowed -\n\
  use `\\' character to override internal functions\n");
      }
      else printff("Error: allow does not take options - operation canceled\n");
    }

    else if(!strcmp(key, "disallow")){
      if(sscanf(str+position, "%s%n", key, &position)!=1){
	ALLOW=0;
	printff("Execution of external (system) commands will not be allowed -\n\
  only internal functions will be recognized\n");
      }
      else printff("Error: disallow does not take options - operation canceled\n");
    }

    else if(!strcmp(key, "susp")){
      saux=str+position;
      if(sscanf(saux, "%s%n", key, &position)!=1){
	suspend=1-suspend;
	printff("Suspended mode: %s\n", onoff[suspend]);
      }

      else if(!strcmp(key, "on")){
	if(sscanf(saux+position, "%s", key)==1)
	  printff("Error: there are more than expected elements in the line\n");
	else{
	  suspend=1;
	  printff("Suspended mode: %s\n", onoff[suspend]);
	}
      }

      else if(!strcmp(key, "off")){
	if(sscanf(saux+position, "%s", key)==1)
	  printff("Error: there are more than expected elements in the line\n");
	else{
	  suspend=0;
	  printff("Suspended mode: %s\n", onoff[suspend]);
	}
      }

      else printff("Error: susp: `%s': unavailable option\n", key);
    }

    else if(!strcmp(key, "exec")){
      if(sscanf(str+position, "%s", key)==1)
	printff("Error: exec: this function takes no parameters\n");
      else{
	printff("Executing now\n");
	if(state==2){
	  if((mass_comp==5)||(mass_comp==6)) fast_flag=1;
	  else fast_flag=0, state=1;
	}
	return 0;
      }
    }

    else if(!strcmp(key, "show")){
      saux=str+position;
      if(sscanf(saux, "%s%n", key, &position)!=1){
	printff("This parser is used to change term number %i for each component\n", mass_comp);
	printff("Number of entries per line is defined to be %i\n\
Element arrangement is shown below:\n", lnum);
	for(i=0;i<lnum;i++) printff("\t%i", mask[i]);
	printff("\nNumber of constant entries is %i\nConstants are:\n", num_mass);
	for(i=0;i<num;i++) printff("\t%g", mass[i]);
	printff("\nSuspended mode is turned %s\n", onoff[suspend]);
	printff("Execution of external (system) commands is %sed\n", allow_disallow[ALLOW]);
	printff("Type `show all' for a list of all variables accessible with this parser\n\n");
      }

      else if(!strcmp(key, "all")){
	if(sscanf(saux+position, "%s", key)==1)
	  printff("Error: only one parameter can be paired with `show'\n");
	else{
	  printff("number of components is %i\n", num);
	  printff("Components are:\n");
	  printff("index real\tindex imag.\tradius avg.\tradius RMS\tconcentration\tdensity\n");
	  for(i=0;i<num;i++){
	    for(j=0;j<NUM_OF_DEF-1;j++) printff("%-15g ", components[i][j]);
	    printff("%g\n", components[i][NUM_OF_DEF-1]);
	  }
	  printff("\n");
	}
      }

      else printff("Error: show: `%s': unavailable option\n", key);
    }

    else if(!strcmp(key, "set")){
      saux=str+position;
      if(sscanf(saux, "%s%n", key, &position)!=1){
	printff("Parameter expected: available options are:\n\
\n\
    `lnum'  sets number of elements per line\n\
    `mask'  sets order of assignment for these elements\n\
    `mass'  sets values of constant terms\n\
    `term'  sets the column from the data file to be changed\n\
\n");
      }

      else if(!strcmp(key, "lnum")){
	saux+=position;
	if(sscanf(saux, "%s%n", key, &position)!=1) printff("Error: No number was provided\n");
	else if(sscanf(saux+position, "%s", key)==1)
	  printff("Error: there are more than expected elements in the line\n");
	else if(sscanf(key, "%i%n", &lnum_aux, &position)!=1)
	  printff("Error: `%s' is not a number\n", key);
	else if(strlen(key)!=position) printff("Error: parsing error\n");
	else if(lnum_aux<num-num_mass) printff("Error: lnum must not be < %i !\n", num-num_mass);
	else if(lnum_aux>MAX_PARSE) printff("Error: lnum must not be > %i !\n", MAX_PARSE);
	else{
	  lnum=lnum_aux;
	  for(i=0, j=0;i<num;i++) if(mass[i]==0){
	    mask[j]=i+1;
	    j++;
	  }
	  for(i=j;i<num;i++) mask[i]=0;
	  for(i=num;i<lnum;i++) mask[i]=0;
	  printff("lnum changed to %i, mask is reset\n", lnum);
	}
      }

      else if(!strcmp(key, "term")){
	saux+=position;
	if(sscanf(saux, "%s%n", key, &position)!=1) printff("Error: No number was provided\n");
	else if(sscanf(saux+position, "%s", key)==1)
	  printff("Error: there are more than expected elements in the line\n");
	else if(sscanf(key, "%i%n", &lnum_aux, &position)!=1)
	  printff("Error: `%s' is not a number\n", key);
	else if(strlen(key)!=position) printff("Error: parsing error\n");
	else if(lnum_aux<=0) printff("Error: term must be > 0 !\n");
	else if(lnum_aux>NUM_OF_DEF) printff("Error: term must not be > %i !\n", NUM_OF_DEF);
	else{
	  mass_comp=lnum_aux;
	  printff("term changed to %i\n", mass_comp);
	}
      }

      else if(!strcmp(key, "mask")){
	saux+=position;
	for(i=0;i<lnum;i++){
	  if(sscanf(saux, "%s%n", key, &position)!=1){
	    printff("Error: there are fewer than expected elements in the line\n");
	    break;
	  }
	  saux+=position;
	  if(sscanf(key, "%i%n", &mask_aux[i], &position)!=1){
	    printff("Error: `%s' is not a number\n", key);
	    break;
	  }
	  if(strlen(key)!=position){
	    printff("Error: parsing error\n");
	    break;
	  }
	}
	if(i==lnum){
	  if(sscanf(saux, "%s", key)==1)
	    printff("Error: there are more than expected elements in the line\n");
	  else{
	    for(j=0;j<=num;j++){
	      flag=0;
	      for(i=0;i<lnum;i++) if(mask_aux[i]==j) flag++;
	      if(j==0) if(flag!=lnum-num+num_mass) break;
	      if(j!=0) if(!mass[j-1]) if(flag!=1) break;
	    }
	    if(j!=num+1) printff("Error: this mask cannot be implemented\n");
	    else{
	      for(i=0;i<lnum;i++) mask[i]=mask_aux[i];
	      printff("Mask changed\n");
	    }
	  }
	}
      }

      else if(!strcmp(key, "mass")){
	saux+=position;
	for(i=0;i<num;i++){
	  if(sscanf(saux, "%s%n", key, &position)!=1){
	    printff("Error: there are fewer than expected elements in the line\n");
	    break;
	  }
	  saux+=position;
	  if(sscanf(key, "%lf%n", &mass_aux[i], &position)!=1){
	    printff("Error: `%s' is not a number\n", key);
	    break;
	  }
	  if(strlen(key)!=position){
	    printff("Error: parsing error\n");
	    break;
	  }
	}
	if(i==num){
	  if(sscanf(saux, "%s", key)==1)
	    printff("Error: there are more than expected elements in the line\n");
	  else{
	    flag=0;
	    for(i=0;i<=num;i++) if(mass_aux[i]!=mass[i]) flag++;
	    if(flag){
	      num_mass=0;
	      for(i=0;i<num;i++){
		mass[i]=mass_aux[i];
		if(mass[i]!=0) num_mass++;
	      }
	      for(i=0, j=0;i<num;i++) if(mass[i]==0){
		mask[j]=i+1;
		j++;
	      }
	      for(i=j;i<num;i++) mask[i]=0;
	      for(i=num;i<lnum;i++) mask[i]=0;
	      printff("Constants are changed, mask is reset\n");
	    }
	    else printff("Constants are kept the same\n");
	  }
	}
      }

      else printff("Error: set: `%s': unavailable option\n", key);
    }

    else if(*key=='#');

    else{
      saux=str;
      position=0;
      for(i=0;i<lnum;i++){
	if(sscanf(saux, "%s%n", key, &position)!=1){
	  printff("Error: there are fewer than expected elements in the line\n");
	  break;
	}
	saux+=position;
	if(sscanf(key, "%lf%n", &mass_aux[i], &position)!=1){
	  if(i==0){
	    if(!ALLOW) printff("Error: keyword `%s' was not found\n", key);
	    else system(str);
	  }
	  else 	  printff("Error: `%s' is not a number\n", key);
	  break;
	}
	if(strlen(key)!=position){
	  printff("Error: parsing error\n");
	  break;
	}
      }
      if(i==lnum){
	if(sscanf(saux, "%s", key)==1)
	  printff("Error: there are more than expected elements in the line\n");
	else{
	  for(i=0;i<lnum;i++) if(mask[i]) components[mask[i]-1][mass_comp-1]=mass_aux[i];
	  for(i=0;i<num;i++) if(mass[i]) components[i][mass_comp-1]=mass[i];

	  /*
	    Enter here any additional modifications you may have,
	    such as linear dependence of components on each other:
	    #define LIN_FACTOR 10.
	    components[3][mass_comp-1]=components[0][mass_comp-1]/LIN_FACTOR;

	    A more versatile approach would be to use mass[i] instead of LIN_FACTOR:
	    components[3][mass_comp-1]=components[0][mass_comp-1]/mass[3];

	    Of course, you can always use awk to take care of these
	    modifications, or to disregard this parser altogether.
	  */

	  if(!NO_MESSAGES){
	    printff("Component values are:\n");
	    for(i=0;i<num;i++) printff("\t%g", components[i][mass_comp-1]);
	    printff("\n");
	  }
	  if(state==2){
	    if((mass_comp==5)||(mass_comp==6)) fast_flag=1;
	    else fast_flag=0, state=1;
	  }
	  if(!suspend) return 0;
	}
      }
    }
  }
}


int coefficient(float *a1, float *b1, float *g1, float *s11, __complex__ float m, float x)
{
  long int i;
  long int k_stop, k_max;
  int error;
  float aux;
  __complex__ float cx_aux;
  __complex__ float mx;
  float cosx, sinx;
  __complex__ float *dzi, dzi_1;
  __complex__ float *D;
  __complex__ float *a, *b;
  __complex__ float *pi, *tu;
  float Q_ext, Q_sca, gG, S11;
  float _Q_ext, _Q_sca, _gG;
  __complex__ float S1, S2, _S1, _S2;

  mx=m*x;
  k_stop=x+4*pow(x, 1./3)+2;
  aux=~mx*mx;
  aux=sqrt(aux);
  if(aux<k_stop) aux=k_stop;
  k_max=aux+15;

  if(DEBUG&128) printff("\nk_max=%i, k_stop=%i, aux=%g\n", k_max, k_stop, aux);

  if(k_max>=K_MAX){
    printf("Unrecoverable error: k_max=%i<K_MAX=%i\n", k_max, K_MAX);
    if(DEBUG&128){
      k_max=K_MAX-1;
      k_stop=k_max-15;
    }
    else return 1;
  }

  error=1;
  if((dzi=malloc((k_max+1)*sizeof(__complex__ float)))==NULL);
  else if((D=malloc((k_max+1)*sizeof(__complex__ float)))==NULL) free(dzi);
  else if((a=malloc((k_max+1)*sizeof(__complex__ float)))==NULL) free(dzi), free(D);
  else if((b=malloc((k_max+1)*sizeof(__complex__ float)))==NULL) free(dzi), free(D), free(a);
  else if((pi=malloc((k_max+1)*sizeof(__complex__ float)))==NULL) free(dzi), free(D), free(a), free(pi);
  else if((tu=malloc((k_max+1)*sizeof(__complex__ float)))==NULL) free(dzi), free(D), free(a), free(pi), free(tu);
  else error=0;
  if(error){
    printf("Unrecoverable error: Memory allocation error\n");
    return 2;
  }

  cosx=cos(x), sinx=sin(x);
  __real__ dzi_1=cosx, __real__ dzi[0]=sinx;
  __imag__ dzi_1=sinx, __imag__ dzi[0]=-cosx;

  dzi[1]=dzi[0]/x-dzi_1;

  for(i=1;i<=k_stop;i++) dzi[i+1]=(2*i+1)*dzi[i]/x-dzi[i-1];

  pi[0]=0, pi[1]=1;
  tu[0]=0; tu[1]=ang;

  for(i=2;i<=k_stop;i++) pi[i]=((2*i-1.)/(i-1))*ang*pi[i-1]-(i/(i-1.))*pi[i-2];
  for(i=2;i<=k_stop;i++) tu[i]=i*ang*pi[i]-(i+1)*pi[i-1];

  if(DEBUG&1){
    printff("\n");
    printff("Step -1: dzi=%g %gi\n", __real__ dzi_1, __imag__ dzi_1);
    for(i=0;i<=k_stop;i++)
      printff("Step %i: dzi=%g %gi\n", i, __real__ dzi[i], __imag__ dzi[i]);
  }

  D[k_max]=0+0i;
  for(i=k_max;i>1;i--) D[i-1]=i/mx-1/(D[i]+i/mx);

  if(DEBUG&2){
    printff("\n");
    for(i=k_max;i>0;i--)
      printff("Step %i: D=%g %gi\n", i, __real__ D[i], __imag__ D[i]);
  }

  for(i=1;i<=k_stop+1;i++){
    cx_aux=D[i]/m+i/x;
    a[i]=(cx_aux*(__real__ dzi[i])-(__real__ dzi[i-1]))/(cx_aux*dzi[i]-dzi[i-1]);
    cx_aux=m*D[i]+i/x;
    b[i]=(cx_aux*(__real__ dzi[i])-(__real__ dzi[i-1]))/(cx_aux*dzi[i]-dzi[i-1]);
  }

  if(DEBUG&4){
    printff("\n");
    for(i=1;i<=k_stop+1;i++)
      printff("Step %i: a=%g %gi, b=%g %gi\n", i,\
	     __real__ a[i], __imag__ a[i], __real__ b[i], __imag__ b[i]);
  }

  if(DEBUG&8) printff("\n");

  Q_ext=0, Q_sca=0, gG=0, S1=0, S2=0;
  for(i=1;i<=k_stop;i++){
    _Q_ext=(2*i+1)*(__real__ a[i] + __real__ b[i]);
    _Q_sca=(2*i+1)*(__real__ (a[i]*~a[i]) + __real__ (b[i]*~b[i]));
    _gG=((i*(i+2.))/(i+1))*(__real__ (a[i]*~a[i+1]) + __real__ (b[i]*~b[i+1]))+\
      ((2*i+1.)/(i*(i+1)))*(__real__ (a[i]*~b[i]));
    _S1=((2*i+1.)/(i*(i+1)))*(a[i]*pi[i]+b[i]*tu[i]);
    _S2=((2*i+1.)/(i*(i+1)))*(a[i]*tu[i]+b[i]*pi[i]);

    if(DEBUG&8)
      printff("Step %i: Q_ext+=%g, Q_sca+=%g, gG+=%g\n", i, _Q_ext, _Q_sca, _gG);

    Q_ext+=_Q_ext, Q_sca+=_Q_sca, gG+=_gG, S1+=_S1, S2+=_S2;
  }

  if(DEBUG&8)
    printff("Sum: Q_ext=%g, Q_sca=%g, gG=%g\n", Q_ext, Q_sca, gG);

  Q_ext=(2/(x*x))*Q_ext;
  Q_sca=(2/(x*x))*Q_sca;
  gG=(4/((x*x)*Q_sca))*gG;
  S11=(S1*~S1+S2*~S2)/(2*(x*x)*Q_sca);

  S11/=Pi; // required for normalization to 1; where does this come from?

  if(DEBUG&16)
    printff("Q_ext=%g, Q_sca=%g, gG=%g\n", Q_ext, Q_sca, gG);

  *a1=Q_ext-Q_sca, *b1=Q_sca, *g1=gG, *s11=S11;

  free(dzi), free(D), free(a), free(b), free(pi), free(tu);
  return 0;
}


/*
  a, b, g:
  absorption and scattering coefficients, asymmetry parameter;
  lambda, n_ice;
  index real, index imag., radius avg., radius RMS.
*/
int integrate(float *a, float *b, float *g, float *s, float lambda, float n_ice,
	      float n_r, float n_i, float r_av, float r_si)
{
  int i;
  int error=0;
  __complex__ float m;
  float x;
  float k_lambda;
  float r;
  float a1, b1, g1, s11;
  float aux, area;

  error=0;
  __real__ m=n_r;
  __imag__ m=n_i;
  m=m/n_ice;
  k_lambda=2*Pi*n_ice/lambda;

  *a=0, *b=0, *g=0, *s=0;

  if(DEBUG&64){
    r=r_av;
    x=k_lambda*r;
    error=coefficient(&a1, &b1, &g1, &s11, m, x);
    printff("\na=%g, b=%g, g=%g, s=%g\n", a1, b1, g1/b1, s11/b1);
    return 0;
  }

  for(i=0;i<=MAX_STEPS;i++){
    r=r_av*pow(r_si, (SIGMA_RANGE_MINUS+SIGMA_RANGE_PLUS)*i/MAX_STEPS-SIGMA_RANGE_MINUS);
    x=k_lambda*r;

    if(DEBUG){
      printff("\n\nStep #%i: proceed?  ", i);
      while(getchar()!='\n');
    }

    area=Pi*(r*r);
    area*=(SIGMA_RANGE_MINUS+SIGMA_RANGE_PLUS)*log(r_si)/MAX_STEPS;
    if((i==0)||(i==MAX_STEPS)) area/=2;
    aux=log(r/r_av)/log(r_si);
    area*=(1/(sqrt(2*Pi)*log(r_si)))*exp(-aux*aux/2);

    error=coefficient(&a1, &b1, &g1, &s11, m, x);

    if(DEBUG&32) printff("\na+=%g, b+=%g, g+=%g, s+=%g\n", a1*area, b1*area, g1*area, s11*area);

    *a+=a1*area;
    *b+=b1*area;
    *g+=b1*area*g1;
    *s+=b1*area*s11;
  }
  *g/=*b;
  *s/=*b;
  return 0;
}


int process(char *filename, float lambda)
{
  int i, j;
  int error=0;
  int num_of_comp=0;
  float n_ice=1., density_ice=1.;
  float components[MAX_NUM_OF_COMP][NUM_OF_DEF];
  float indexf[5];
  float a, b, g, s, at, bt, gt, st;
  float abg[MAX_NUM_OF_COMP][4];
  float aux;
  char str[MAX_LENGTH_FLINE];
  char ch;

  FILE *fp;
  if((fp=fopen(filename,"r"))==NULL){
    printff("cannot open output file `%s', exiting\n", filename);
    return 1;
  }
  error=0;

  if(fscanf(fp, "%lf", &n_ice)!=1) error=1;
  if(fscanf(fp, "%lf", &density_ice)!=1) error=1;
  if(fscanf(fp, "%i", &num_of_comp)!=1) error=1;
  if(!error) if(num_of_comp<0)
    printff("Nothing to be done for `%s'\n", filename), error=1;
  if(!error) if(num_of_comp>MAX_NUM_OF_COMP)
    printff("Maximum allowed number of components exceeded\n"), error=1;
  if(!error) for(i=0;i<num_of_comp;i++) for(j=0;j<NUM_OF_DEF;j++){
    if(fscanf(fp, "%lf", &components[i][j])!=1) error=1;
    if(error) break;
  }
  if(USE_INDEX) if(!error){
    if(!NO_MESSAGES) printff("\nAdditional information for refractive index reconstruction:\n\
c0\tc1\tc2\td0\tb\n");
    aux=log(lambda);
    if(!NO_MESSAGES) printff("Ice:\n");
    for(j=0;j<3;j++){
      if(fscanf(fp, "%lf", &indexf[j])!=1) error=1;
      if(error) break;
      if(!NO_MESSAGES) printff("%g\t", indexf[j]);
    }
    n_ice=indexf[0]+indexf[1]*aux+indexf[2]*aux*aux;

    if(!NO_MESSAGES) printff("\nImpurity components:\n");
    for(i=0;i<num_of_comp;i++){
      for(j=0;j<5;j++){
	if(fscanf(fp, "%lf", &indexf[j])!=1) error=1;
	if(error) break;
	if(!NO_MESSAGES) printff("%g\t", indexf[j]);
      }
      if(error) break;
      components[i][0]=indexf[0]+indexf[1]*aux+indexf[2]*aux*aux;
      components[i][1]=indexf[3]*pow((lambda/337), indexf[4]);
      if(!NO_MESSAGES) printff("\n");
    }
  }
  if(error){
    printff("file format error\n");
    fclose(fp);
    return 1;
  }

  printff("refractive index of ice is %g\n", n_ice);
  printff("density of ice is %g\n", density_ice);
  printff("wavelength of the incident light is %g nm\n", lambda);
  printff("number of components is %i\n", num_of_comp);
  printff("index real\tindex imag.\tradius avg.\tradius RMS\tconcentration\tdensity\n");
  for(i=0;i<num_of_comp;i++){
    for(j=0;j<NUM_OF_DEF-1;j++) printff("%-15g ", components[i][j]);
    printff("%g\n", components[i][NUM_OF_DEF-1]);
  }

  do{
    if(fscanf(fp, "%c", &ch)!=1) error=1;
  } while((ch!='#')&&(!error));
  do{
    error=0;
    for(i=0;i<MAX_LENGTH_FLINE-1;i++){
      if(fscanf(fp, "%c", &ch)!=1) error=1;
      if(ch!='\n') str[i]=ch;
      else break;
    }
    if(i==MAX_LENGTH_FLINE-1) do{
      if(fscanf(fp, "%c", &ch)!=1) error=1;
    } while((ch!='\n')&&(!error));
    str[i]='\0';

    if(!NO_MESSAGES) printff("%s\n", str);
  } while(!error);

  fclose(fp);

  while(!input(components, num_of_comp)){
    at=0, bt=0, gt=0, st=0;
    if(!NO_MESSAGES){
      printff("\
Absorption and scattering coefficients, and asymmetry parameter found are:\n\
a, [m^-1]\tb, [m^-1]\tg = <cos(O)>\ts11\t\t");
      if(INCLUDE_ICE) printff("b_e=b (1-g)\t");
      printff("n, [nm^-3]\n");
    }
    for(i=0;i<num_of_comp;i++){

      if(!fast_flag){
	error=integrate(&a, &b, &g, &s, lambda, n_ice, components[i][0],
			components[i][1], components[i][2]*1000, components[i][3]);
	if(error) return 1;
	abg[i][0]=a, abg[i][1]=b, abg[i][2]=g, abg[i][3]=s;
      }
      else a=abg[i][0], b=abg[i][1], g=abg[i][2], s=abg[i][3];

      aux=components[i][4]/((4*Pi/3)*(components[i][5]/density_ice)*\
			    (components[i][2]*components[i][2]*components[i][2]*1.e9)*\
			    exp(4.5*log(components[i][3])*log(components[i][3])));
      a*=aux, b*=aux;
      if(!NO_MESSAGES){
	printff("%-15g %-15g %-15g %-15g ", a, b, g, s);
	if(INCLUDE_ICE) printff("%-15g ", b*(1-g));
	printff("%g\n", aux);
      }
      at+=a, bt+=b, gt+=g*b, st+=s*b;
    }
    gt/=bt, st/=bt;
    if(!NO_MESSAGES){
      printff("Final summed values:\n");
      if(INCLUDE_ICE) printff("a_all\t\tb_e\t\t");
      printff("a\t\tb\t\tg\t\ts11\n");
    }
    if(INCLUDE_ICE) printf("%-15g %-15g ", at+ICE_A*exp(-ICE_al*lambda)+ICE_B*exp(-ICE_bl/lambda), bt*(1-gt));
    printf("%-15g %-15g %-15g %-15g\n", at, bt, gt, st);
    if(!VERBAL) break;
  }

  return 0;
}


int main(int arg_c, char *arg_a[])
{
  static int i=0, flag=0;
  float lambda=LAMBDA_DEFAULT, cth=10;
  float SR_M=SIGMA_RANGE_MINUS, SR_P=SIGMA_RANGE_PLUS;
  float STEPS=MAX_STEPS;

  char *ex_body = "\
1.31 0.917 4\n\
\n\
1.53 0.005 .27 2.67 12 2.6\n\
1.39 .3e-8 .4 2.03 33 2.2\n\
1.43 1.e-8 .0695 1.86 171 1.68\n\
1.75 0.4 .0118 2.00 .4 2.3\n\
\n\
2.544 -0.3699 0.02764\n\
\n\
1.530 0 0 0.0188 -2.448\n\
2.692 -0.4117 0.0323 1.3e-7 -8.278\n\
6.112 -1.4907 0.1187 1.0e-8 0\n\
-0.516 0.7418 -0.0607 0.467 -0.102\n\
\n\
***comments***\n\
#order of components:\n\
mineral salt acids soot\n\
\n";

  char *description = "\
This program computes scattering coefficients using Mie theory, v%g\n\
Available options are: -h  = -help  print this message and quit\n\
                       -q  = -quit  quit immediately\n\
                       -L=[wavelength of incident light in nm]\n\
                       -SR-=[SIGMA_RANGE_MINUS]  range of integration\n\
                             below average radius, in ln(GSD)'s\n\
                       -SR+=[SIGMA_RANGE_PLUS]   same above it\n\
                             default is SR-= %g, SR+= %g\n\
                       -STEPS=[MAX_STEPS] number of integration steps\n\
                             default is STEPS=%i\n\
                       -n    use formulae for the refractive index\n\
                            (additional info in the input file needed)\n\
                       -u = -usage description of data sets and output\n\
                       -d[debug code] - invokes step by step execution\n\
                       -sm   suppress auxiliary messages\n\
                       -zm   kill all messages, but the result line\n\
                       -v    ask for mass densities of the components\n\
                       -ice  show a_all, that includes a contribution\n\
                             to a(L) from the ice, and b_e=b(1-g)\n\
                       -th=[scattering angle in degrees]  for S11\n\
                       -ct=[cos of the scattering angle]  for S11\n\
                       -ex   print an example of an input file\n\
                       -a    allow execution of external commands\n\
\n";
  char *usage ="\
Data file has the following structure:\n\
First 3 numbers are: refractive index and density of ice, and number of\n\
components to be taken from the list to follow. Each line in this list has\n\
%i elements: real and imaginary parts of the refractive index, mean radius\n\
(in mkm), and its geometrical standard deviation (GSD), mass concentration\n\
(in ng/g), and density of the component. It is assumed that particle radii\n\
are distributed log-normally with the mean radius and GSD defined above.\n\
\n\
All densities must be expressed in the same units, e.g. [g/cm^3].\n\
Wavelength of the incident light must be expressed in nm via option -L,\n\
by default it is equal to %g nm.\n\
\n\
On the output absorption (a) and scattering (b) coefficients are expressed\n\
in 1/m, and asymmetry parameter (g=<cos theta>) is dimensionless.\n\
\n\
If option -n is specified, additional section in the data file is needed.\n\
For each component a line, containing 5 elements c0, c1, c2, d0, b will be\n\
read in. First line is for the ice, so only 3 components are needed, since\n\
refractive index of the ice is assumed to be purely real.\n\
In such a case information about the refractive indices, contained in\n\
the previous section of the file, will be disregarded, and the following\n\
formula will be used instead: N = c0 + c1 ln L + c2 ln^2 L + i d0 (L/337)^b\n\
\n\
Option -ice shows a_all = %g exp-%gL + %g exp-%g/L + a(L)\n\
and b_e = b(L) [1-g(L)]. First two terms in the expression for a_all\n\
describe contribution of the pure ice to the absorption.\n\
\n";

  for(i=1;i<arg_c;i++){
    if(*arg_a[i]=='-'){

      if(!strcmp(arg_a[i], "-ice")) INCLUDE_ICE=1;
      else if(!strcmp(arg_a[i], "-v")) VERBAL=1;
      else if(!strcmp(arg_a[i], "-sm")) NO_MESSAGES=1;
      else if(!strcmp(arg_a[i], "-zm")) ZERO_MESSAGES=1;
      else if(!strcmp(arg_a[i], "-n")) USE_INDEX=1;
      else if(!strcmp(arg_a[i], "-a")) ALLOW=1;
      else if(!strncmp(arg_a[i], "-STEPS=", 7)) STEPS=atof(arg_a[i]+7);
      else if(!strncmp(arg_a[i], "-SR-=", 5)) SR_M=atof(arg_a[i]+5);
      else if(!strncmp(arg_a[i], "-SR+=", 5)) SR_P=atof(arg_a[i]+5);
      else if(!strncmp(arg_a[i], "-L=", 3)) lambda=atof(arg_a[i]+3);
      else if(!strncmp(arg_a[i], "-th=", 4)) ang=atof(arg_a[i]+4);
      else if(!strncmp(arg_a[i], "-ct=", 4)) cth=atof(arg_a[i]+4);
      else if(!strncmp(arg_a[i], "-d", 2)) DEBUG=atoi(arg_a[i]+2);

      else if(!strcmp(arg_a[i], "-ex")){
	printf("%s", ex_body);
	return 1;
      }
      else if(!(strcmp(arg_a[i], "-u")&&strcmp(arg_a[i], "-usage"))){
	printf(usage, NUM_OF_DEF, LAMBDA_DEFAULT, ICE_A, ICE_al, ICE_B, ICE_bl);
	return 1;
      }
      else if(!(strcmp(arg_a[i], "-q")&&strcmp(arg_a[i], "-quit"))) return 1;
      else if(!(strcmp(arg_a[i], "-h")&&strcmp(arg_a[i], "-help"))){
	printf(description, MYVERSION, SIGMA_RANGE_MINUS, SIGMA_RANGE_PLUS, MAX_STEPS);
	return 1;
      }
      else{
	printf("Unknown option: `%s', aborting execution\n", arg_a[i]);
	return 1;
      }
    }
    else flag=1;
  }

  if(flag){
    if(SR_P>-SR_M) SIGMA_RANGE_MINUS=SR_M, SIGMA_RANGE_PLUS=SR_P;
    else printff("Warning: The lower limit of integration is greater than the \
higher,\nlimits reversed to defaults\n");
    if(STEPS>0) MAX_STEPS=STEPS;
    else printff("Warning: Number of steps zero or negative, reversing to default\n");
    if(lambda<=0) lambda=LAMBDA_DEFAULT, printff("Warning: L must be > 0, L is reset to default value\n");
    if(cth>=-1 && cth<=1) ang=cth; else ang=cos(ang*Pi/180);
    if(NO_MESSAGES) printff("Auxiliary messages will be suppressed\n");
    else{
      if(DEBUG) printff("Debug option %i enabled\n", DEBUG);
      printff("number of integration steps is %i\n", MAX_STEPS);
      printff("integration goes from -%g*sigma to %g*sigma\n",
	     SIGMA_RANGE_MINUS, SIGMA_RANGE_PLUS);
      if(INCLUDE_ICE) printff("a_all and b_e will be shown\n");
      else printff("a_all and b_e will not be shown: use option -ice\n");
      if(USE_INDEX) printff("Index formulae will be used,\n\
additional information in the data file is required\n");
      if(VERBAL){
	printff("You will need to supply mass densities at the program prompt\n");
	if(ALLOW) printff("Use of external commands will be recognized\n");
      }
    }

    for(i=1;i<arg_c;i++) if(*arg_a[i]!='-'){
      printff("file `%s' is being fed in\n", arg_a[i]);
      if(process(arg_a[i], lambda)) printff("error in function `process'\n");
    }
  }

  else{
    printff("Usage: %s [option(s)] [input file(s)]\n\
Type %s -h for help and a list of available options\n", arg_a[0], arg_a[0]);
    return 1;
  }

  return 0;
}