#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include "complex.c"
#define Pi 3.14159265
#define NANG  90
#define REFMED 1.0
#define NWAVE 21
/* parameters for an ice cloud */

double *bhmie(double,fcomplex,int);
double compute_qsca(double);
double compute_qabs(double);
double compute_g(double);

float lambda,refre,refim;

main(int argc,char **argv)
{
  FILE *indexfp;
  double rad;
  int istring,iwave;
    for(rad = 1.;rad < 20. ; rad += 9.){
    for(lambda = .5;lambda <= 20.;lambda += 1000.) {
    refre = 1.44;refim = 0.;
    printf("%f\t%f\t%f\t%f\t%f\t",rad,lambda,compute_qsca(rad),
    compute_qabs(rad),compute_g(rad));
  }printf("\n");}
}

double compute_qsca(double rad)
{
  fcomplex refrel;
  double x,*q,qsca;

  refrel = RCmul(1/REFMED,Complex(refre,refim));
  x = 2*Pi*rad*REFMED/lambda;
  q = bhmie(x,refrel,NANG);
  qsca = *(q + 1);

  return qsca;
}

double compute_qabs(double rad)
{
  fcomplex refrel;
  double x,*q,qext,qsca,qabs;

  refrel = RCmul(1/REFMED,Complex(refre,refim));
  x = 2*Pi*rad*REFMED/lambda;
  q = bhmie(x,refrel,NANG);
  qext = *q;
  qsca = *(q + 1);
  qabs = qext - qsca;

  return qabs;
}

double compute_g(double rad)
{
  fcomplex refrel;
  double x,*q,g;

  refrel = RCmul(1/REFMED,Complex(refre,refim));
  x = 2*Pi*rad*REFMED/lambda;
  q = bhmie(x,refrel,NANG);
  g = *(q + 2);

  return g;
}

double *bhmie(double x,fcomplex refrel,int nang)
{
  double amu[100],theta[100],pi[100],tau[100],pi0[100],pi1[100],
         psi0,psi1,psi,dn,dx;
  fcomplex d[3000],y,xi,xi0,xi1,an,bn,an_old,bn_old,s1[200],s2[200];  
  double xstop,dang,ymod,chi0,chi1,apsi0,apsi1,rn,fn,apsi,chi;
  int nstop,nmx,j,n,nn,jj,p,t;
  double qext,qsca,qback,q[3],g;
  
  dx = x;
  y = RCmul(x,refrel);

  xstop = x + 4*pow(x,1/3) + 2;
  nstop = xstop;
  ymod = Cabs(y);
  nmx = (xstop > ymod) ? xstop + 15 : ymod + 15;
  dang = Pi/2/(double)(nang - 1);

  for(j = 1;j <= nang;j++)
  {
    theta[j] = ((double)j - 1)*dang;
    amu[j] = cos(theta[j]);
  }

  d[nmx] = Complex(0.0,0.0);
  nn = nmx - 1;
  for(n = 1;n <= nn;n++)
  {
    rn = nmx - n + 1;
    d[nmx-n] = Csub(Cdiv(Complex(rn,0),y),
		    Cdiv(Complex(1,0),Cadd(d[nmx-n+1],Cdiv(Complex(rn,0),y))));
  }

  for(j = 1;j <= nang;j++)
  {
    pi0[j] = 0.0;
    pi1[j] = 1.0;
  }
  nn = 2*nang - 1;
  for(j = 1;j <= nn;j++)
  {
    s1[j] = Complex(0.0,0.0);
    s2[j] = Complex(0.0,0.0);
  }

  psi0 = cos(dx);
  psi1 = sin(dx);
  chi0 = -sin(x);
  chi1 = cos(x);
  apsi0 = psi0;
  apsi1 = psi1;
  xi0 = Complex(apsi0,-chi0);
  xi1 = Complex(apsi1,-chi1);
  qsca = 0.0;
  g = 0.0;
  n = 1;
  do
  {
    dn = n;
    rn = n;
    fn = (2*rn + 1)/(rn*(rn + 1));
    psi = (2*dn - 1)*psi1/dx - psi0;
    apsi = psi;
    chi = (2*rn - 1)*chi1/x - chi0;
    xi = Complex(apsi,-chi);
    an = Csub(RCmul(apsi,Cadd(Cdiv(d[n],refrel),Complex(rn/x,0))),
	      Complex(apsi1,0));
    an = Cdiv(an,Csub(Cmul(Cadd(Cdiv(d[n],refrel),Complex(rn/x,0)),xi),xi1));
    bn = Csub(RCmul(apsi,Cadd(Cmul(refrel,d[n]),Complex(rn/x,0))),
	      Complex(apsi1,0));
    bn = Cdiv(bn,Csub(Cmul(Cadd(Cmul(refrel,d[n]),Complex(rn/x,0)),xi),xi1));
    qsca+=(2*rn + 1)*(pow(Cabs(an),2) + pow(Cabs(bn),2));
    if(rn > 1)
      g+=(rn - 1)*(rn + 1)/rn*
 	 (Cadd(Cmul(an_old,Conjg(an)),Cmul(bn_old,Conjg(bn)))).r +
	 (2*(rn - 1) + 1)/(rn - 1)/rn*(Cmul(an_old,Conjg(bn_old))).r;
    an_old = an;
    bn_old = bn;

    for(j = 1;j <= nang;j++)
    {
      jj = 2*nang - j;
      pi[j] = pi1[j];
      tau[j] = rn*amu[j]*pi[j] - (rn + 1)*pi0[j];
      p = pow(-1,n - 1);
      s1[j] = Cadd(s1[j],RCmul(fn,Cadd(RCmul(pi[j],an),RCmul(tau[j],bn))));
      t = pow(-1,n);
      s2[j] = Cadd(s2[j],RCmul(fn,Cadd(RCmul(tau[j],an),RCmul(pi[j],bn))));
      if(j == jj) continue;
      s1[jj] = Cadd(s1[jj],
		    RCmul(fn,Cadd(RCmul(pi[j]*p,an),RCmul(tau[j]*t,bn))));
      s2[jj] = Cadd(s2[jj],
		    RCmul(fn,Cadd(RCmul(tau[j]*t,an),RCmul(pi[j]*p,bn))));
    }
    
    psi0 = psi1;
    psi1 = psi;
    apsi1 = psi1;
    chi0 = chi1;
    chi1 = chi;
    xi1 = Complex(apsi1,-chi1);
    n = n + 1;
    rn = n;

    for(j = 1;j <= nang;j++)
    {
      pi1[j] = (2*rn - 1)/(rn - 1)*amu[j]*pi[j];
      pi1[j] = pi1[j] - rn*pi0[j]/(rn - 1);
      pi0[j] = pi[j];
    }
  }
  while(n - 1 - nstop < 0);

  qsca*=2/pow(x,2);
  qext = 4/pow(x,2)*s1[1].r;
  qback = 4/pow(x,2)*pow(Cabs(s1[2*nang - 1]),2);
  g*=4/pow(x,2)/qsca;

  q[0] = qext;
  q[1] = qsca;
  q[2] = g;
  return q;
}

#undef NPOINTQUAD
#undef Pi
#undef REFMED
#undef NWAVE
#undef NANG
#undef RMIN
#undef RMAX
#undef N
#undef MODERAD
#undef SIGMA