c gd.f: driver program for g.f
c
c 13 jan: separated subroutines into gsub.f, to work on driver
c
c main program.(version in which eigenvalues are found from matrix)
      program spheroid
      implicit real*8 (a-h,o-z)
      common /accura/ epslon,epslo2
      complex*16 cfindx
      common/param/cfindx,rfindx,anglin,y,lsig
      complex*16 cin
      common/cincou/ cin,cout,sgn,zeta
      common/lwlim/ asize,bsize,wave 
      logical needed
      logical wanted
      common/output/ wanted
      common/conver/ cnvrgn, cnvrgs
      complex*16 dcmplx
      real*8 jbs,nbs
      complex*16 jbsi,nbsi
      common/spfnct/p(120,22),q(120,22),pl(120),jbs(120),nbs(120),
     $     jbsi(120),nbsi(120),qneg(22),maxdeg,maxord
      character*24 sfnamec,sfname
c     
c some hardwired common input and constants
      data epslon/.5d-15/
c
      epslo2=epslon*epslon
      pi=3.141592653589793d0
      maxord=22
      maxdeg=120
c
c read common input
c
      read(5,*) wanted
      read(5,*) needed
      read(5,*) cnvrgn, cnvrgs
      write(6,879) cnvrgn, cnvrgs
 879  format(1x,' convergence criteria = ',2(1pd11.3,2x))
      read(5,*) dangli
      anglin=pi/180.0d0*dangli
      read(5,*) sgn
      read(5,*) aoverb,lsig,y
c lsig=1 if zeta is near 1.0 in prolate case   input also y=zeta-1
c
      zeta=1.0d0/dsqrt(1.0d0-1.0d0/aoverb/aoverb)
      if(sgn.eq.-1.0d0) zeta=dsqrt(1.0d0/(1.0d0-1.0d0/aoverb/aoverb)-1.
     $     0d0)
c
c whenever zeta is updated, need new functions plegen, qlegen, qhyper
      mold=0
c
      if(sgn.eq.1.0d0) then
c prolate options
         write(6,779) dangli,aoverb
 779     format(1x,'prolate static with inclination',f7.3,
     $        ' degrees, a/b = ',f7.4)
      else
         write(6,777) dangli,aoverb
 777     format(1x,'oblate inclined by',f7.3,' degrees, a/b = ',f7.4)
      endif
c
c end of common input
c
c loops:  size, wavelength, trials: change order (and dataset) as appropriate
      read(5,*) wave,rfindx,findxi
      read(5,*) nsize,sizebot,sizedif
      read(5,*) ntry,kmaxpup,mtopup
c
c open summary files
      iaoverb=idint(10.*aoverb)
      irfindx=idint(1000.*rfindx)
      if(sgn.eq.1.d0) then
         write(sfnamec,1110) iaoverb,irfindx
 1110    format('sx',i2,'p',i4,'out.datcheck')
         write(sfname,1111) iaoverb,irfindx
 1111    format('sx',i2,'p',i4,'out.dat')
         write(6,1111) iaoverb,irfindx
      else
         write(sfnamec,1108) iaoverb,irfindx
 1108    format('sx',i2,'o',i4,'out.datcheck')
         write(sfname,1109) iaoverb,irfindx
 1109    format('sx',i2,'o',i4,'out.dat')
         write(6,1109) iaoverb,irfindx
      endif
      open(unit=16,file=sfnamec,status='new')
      open(unit=17,file=sfname,status='new')
      write(16,500)
      write(17,500)
 500  format(3x,'bk',9x,'area*k**2',2x, 'qexte',6x, 'qextm',6x,
     $     'tau',8x,'polar',6x, 'qpte',7x, 'qptm',7x,'birefr',5x,
     $     'qscte',6x, 'qsctm')
c
c loop size ::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
      do 11 isize=1,nsize
         bsize=sizebot+sizedif*dfloat(isize-1)
c wave and size in same units (usually microns)
c bsize is always the short dimension
      asize=aoverb*bsize
      write(6,*) ' '
      write(6,607) asize,bsize
 607  format(1x,'asize = ',f7.4,' bsize =',f7.4)
c end of specific size input :::::::::::::::::::::::::::::::::::::::::
c
c loop trials ****************************************************
      do 1 itry=1,ntry
      write(6,*) ' '
      if(findxi.eq.0.d0) write(6,650)
 650  format(1x,'warning - real index - consider other program')
cm* m = n + ik
cm* note positive sign of - opposite to van de hulst, like born and wolf
      findxi = dabs(findxi)
      cfindx =dcmplx(rfindx,findxi)
c end of specific wavelength information ******************************
c
      scaleb=bsize/wave
      scalea=asize/wave
      ak=2.0d0*pi*scalea
      bk=2.0d0*pi*scaleb
      rhob=2.0d0*(rfindx-1.0d0)*bk
c
      cout=bk*dsqrt(aoverb**2-1.d0)
c  wronskian for outside c
      w=1.0d0/cout/ (zeta*zeta-sgn)
      if(wanted) write(6,606) zeta,w
 606  format(1x,5hzeta=,d23.16,2x,10hwronskian=,d23.16)
      cin=cout*cfindx
c
      kmaxp=kmaxpup - 2*(itry-1)
      mbot=1
      mtop=mtopup
c 
c pair of c,kmaxp-mtop cards with largest mtop should be placed first
c
c note ::::****::::**** ::::****::::**** ::::****::::**** ::::****::::**** 
c consider a new choice every time wavelength and/or size changes
c note ::::****::::**** ::::****::::**** ::::****::::**** ::::****::::**** 
      if(anglin.eq.0.0d0) mbot=2
      if(anglin.eq.0.0d0) mtop=2
      mmax=mtop+1
      nmax=mtop+99
      if(mmax.gt.22) then
         write(6,*) ' maximum value of mtop is 21 :',mtop
         stop
      endif
c
c maximum value for mtop is 21 ; maxord=22
c 21+99=120 which is the other magic dimension, maxdeg
c
      if(kmaxp.gt.39) then
         write(6,*) ' maximum value of kmaxp is 39 :',kmaxp
         stop
      endif
c
      if(mod(kmaxp,2).eq.0) then
         write(6,*) ' kmaxp should be odd :',kmaxp
         stop
      endif
c
      write(6,608) wave,cfindx,kmaxp,mtop
 608  format(1x,'wave =',f7.4,2x,17hrefractive index=,f6.3,f7.3,
     $     2x,6hkmaxp=,i2,2x,5hmtop=,i2)
      write(6,605) cout,cin,scalea,scaleb,ak,bk,rhob
 605  format(1x,2hc=,f7.4,2x,'cin=',2x,f7.4,f8.4,2x,
     $     2ha=,f7.4,2h*w,2x,2hb=,f7.4,2h*w,
     $     2x,4ha*k=,f7.4,2x,4hb*k=,f7.4,2x,12h2*k*b*(n-1)=,f7.4)
c
c end of reading input
c
c calculate single particle
      call sinsph(needed,mold,nmax,mmax,mbot,mtop,kmaxp,cnvrgn)
c
c output single particle results
c
c check for convergence with respect to kmaxp
c
      if(itry.eq.1) then
         call store(1)
         call sinsum(aoverb,bk,17)
         call sinfin(aoverb)
         call sinsum(aoverb,bk,16)
c
c long wavelength limit approximation if appropriate
         akbk=dmin1(ak,bk)
         if(akbk.lt.5.d-2) call lwl
      else
         call store(2)
         call contry(iftry)
         call sinfin(aoverb)
         call sinsum(aoverb,bk,16)
         if(iftry.eq.0) then
            write(6,102) kmaxpup
 102        format(1x,'converged for kmaxp = ',i2)
          go to 101 
          else
             if(itry.eq.ntry) then
                write(6,103) kmaxpup
 103            format(1x,'failed to converge for kmaxp = ',i2)
             endif 
          endif
       endif
c end of trial loop
 1    continue
c skip from trial loop
 101  continue
c end of size loop
 11   continue
      close(16)
      close(17)
      stop
      end
c
      subroutine sinsum(aoverb,bk,iunit)
c write output from single particle calculation
      implicit real*8 (a-h,o-z)
      complex*16 cexte,cextm,cscate,cscatm
      common/outpu/cexte,cextm,cscate,cscatm,iconv
      complex*16 cfindx
      common/param/cfindx,rfindx,anglin,y,lsig
      complex*16 cin
      common/cincou/ cin,cout,sgn,zeta
c
      pi=3.141592653589793d0
      cee=cexte
      cme=cextm
      cep=(0.0d0,-1.0d0)*cexte
      cmp=(0.0d0,-1.0d0)*cextm
c
c basic information before normalization
c
      sumext=(cextm+cexte)/2.d0
      pol=(cextm-cexte)/2.d0
      phasdf=dcmplx(0.0d0,-1.0d0)*(cextm-cexte)/2.d0
c
c scaled areas (have a k**2 factor, as do the cross-sections)
      ab=aoverb**2
      if(sgn.eq.-1.0d0) ab=1.0d0/ab
      a=sgn*pi*cout*cout/(ab-1.0d0)
      cos2=dcos(anglin)**2
      sin2=dsin(anglin)**2
      proj1=cos2+sin2/ab
      proj2=(sin2*cos2*(1.0d0/ab-1)**2)/(sin2+cos2/ab)
      pa=a/dsqrt(proj1)
      pga=a/dsqrt(proj1-proj2)
      pab=pi*cout*cout*aoverb/(aoverb**2-1.0d0)
c
c scaled volume
c need to input approariate parameters if writing this
c      if(sgn.eq.1.0d0) then
c          volume=4.0d0*pi*scalea*scaleb**2/3.0d0
c          else
c             volume=4.0d0*pi*scaleb*scalea**2/3.0d0
c          endif
c       write(6,505) volume
c505   format(' ',7hvolume=,d10.4,5h*w**3)
c
c choose edge-on (scaled) area
c
      qtepab=cexte/pab
      qtmpab=cextm/pab
      qpepab=dimag(cexte)/pab
      qpmpab=dimag(cextm)/pab
      qscate=cscate/pab
      qscatm=cscatm/pab
      extpab=sumext/pab
      polpab=pol/pab
      birpab=phasdf/pab
      brdcab=birpab*polpab
      write(iunit,504)bk,pab,qtepab,qtmpab,extpab,polpab,
     $     qpepab,qpmpab,birpab,qscate,qscatm
 504  format(1x,11(1x,1pd10.3))
c     
      return
      end
c
      subroutine store(ic)
      implicit real*8 (a-h,o-z)
      complex*16 cexte,cextm,cscate,cscatm
      common/outpu/cexte,cextm,cscate,cscatm,iconv
      common/stor/ cstore(6,2)
      cstore(1,ic)=dreal(cexte)
      cstore(2,ic)=dimag(cexte)
      cstore(3,ic)=dreal(cextm)
      cstore(4,ic)=dimag(cextm)
      cstore(5,ic)=dreal(cscate)
      cstore(6,ic)=dreal(cscatm)
      return
      end
c
      subroutine contry(iftry)
      implicit real*8 (a-h,o-z)
      common/conver/ cnvrgn, cnvrgs
      common/stor/ cstore(6,2)
      iftry=0
      do 1 i=1,6
         diff=2.d0*(cstore(i,1)-cstore(i,2))/(cstore(i,1)+cstore(i,2))
         if(dabs(diff).gt.cnvrgn) iftry=1
 1    continue
      return
      end
c



