c     Program to calculate gyrosynchrotron radiation from a power law
c     electron spectrum with isotropic pitch angle distribution. The
c     program uses the theory developed in Ramaty (ApJ 158, 753, 1969).
c     That paper and a paper (Ramaty et al. 1994 ApJ 436, 941) should be
c     referenced in connection with this program.

c     The program requires an input file, gs.inp; the parameters in
c     this file define: 

c     electron spectrum, N(E) = anor E^(-ex)  (MeV^-1)
c     bmag, the magnetic field in Gauss assumed uniform
c     cs, the cosine of the angle between B and the line of sight (0<cs<1)
c        cs can take on any positive value greater than 0 and less
c        than 1; the appropriate formulae for cs=0 are given in Ramaty 
c        (1969); these formulae have not yet been incorporated into the
c        present routine; it appears that the routine works
c        well for 0.1<cs<0.95.
c     theta angular size of source in arc sec; source assumed spherical
c     j1 and j2 define the range of electron energies used in the 
c        integration, for example
c        j1=1; j2=40 corresponds to 0.01 - 100 MeV
c        j1=1; j2=30 corresponds to 0.01 - 10 MeV
c     kf defines the highest frequency (in units of gyrofrequency), examples
c        kf=101 corresponds to ffb=10
c        kf=111 corresponds to ffb=100
c        kf=121 corresponds to ffb=1000
c     etr (in MeV) is an assumed transition energy;
c        below etr the full gyrosynchrotron formalism is used; above etr 
c        the program uses the ultrarelativistic (Ginzburg and Syrovatskii) 
c        formula; to calculate the total emissions it is sufficient
c        to take etr about 3 MeV, except for cs close to 1 when etr should 
c        be larger        
c     alpha is the Razin parameter (Ramaty 1969) = 
c     1.5 gyrofrequency/plasma frequency

c     the given parameters of the current input file will reproduce
c     the middle curve of figure 12 in Ramaty et al. (1994, ApJ, 436, 941)
c     when column 4 is plotted against column 1, lines 131 to 251 in
c     gs.out

c     e1d,e2d ord/xord mode emissivities in erg s^-1 sr^-1 Hz^-1 cm^-3
c     a1d,a2d ord/xord mode abs coefficients in cm^-1
c        both the emissivities and absorption coefficients are 
c        normalized to an electron number density, E^(-ex) (MeV^-1 cm^-3). 
c     phi1,phi1 ord/xord mode flux densities at Earth in 
c        erg s^-1 Hz^-1 cm^-2 normalized to the total electron number
c        N(E) = anor E^(-ex)  (MeV^-1)
c     rccd and rcd are degrees of polarization calculated 
c     from the emissivities and the flux densities, respectively;
c     All 8 quantities are arrays of size defined by kf.
c
c     As the programs runs, it prints on the screen the current ffb 
c     (frequency over gyrofrequency) at which the calculation is carried 
c     out. The program is very fast at low ffb's; then, depending on
c     ETR it will slow down. To reach ffb=1000 with ETR=10 MeV and
c     j2=30 about 60 seconds are required on a Pentium/450MHz. However, if
c     ETR is lower, the program is much faster.  With this ETR the result 
c     is quite accurate, except the polarization.

c     The output is written into gs.out.
c============================================================
      dimension e1d(200),e2d(200),a1d(200),a2d(200),rccd(200)
      dimension ffbd(200),phi1(200),phi2(200),phit(200),rcd(200)
      dimension q(200),v(200),scd(200)
      open(unit=20,file='gs.inp',status='old')
      open(unit=10,file='gs.out',status='unknown')
c----------------read input----------------------------------
      read(20,200)ex
      read(20,200)anor
      read(20,200)bmag
      read(20,200)cs
      read(20,200)theta
      read(20,201)j1
      read(20,201)j2
      read(20,201)kf
      read(20,200)etr
      read(20,200)alpha
200   format(10x,e15.3)
201   format(10x,i10)
c------------------------------------------------------------
      ss=sqrt(1.-cs*cs)
      th=theta/57.3/60./60.
      ffp=1.5/alpha
c------------start loop on ffb-------------------------------
      do k=1,kf
      if(k.le.100)ffb=10.**(0.01*(k-1))
      if(k.gt.100)ffb=10.**(0.1*(k-91))
      ffbd(k)=ffb
      call refr(ffb,ffp,cs,an1,an2,ath1,ath2)
      write(11,103)ffb,an1,an2,ath1,ath2
103   format(5e13.3)
      e1=0.
      e2=0.
      a1=0.
      a2=0.
c----------start integration over energy--------------------
      do j=j1,j2
      el=10.**(0.1*(j-1)-2)
      eu=10.**(0.1*j-2)
      em=10.**(0.1*(j-0.5)-2)
      de=eu-el
      gamma=em/0.511+1.
      gtr=etr/0.511+1.
      g1=0.
      g2=0.
      if(gamma.lt.gtr)then
        if(ffb.gt.ffp)then
        an=sqrt(an1)
        call gsy(gamma,ffb,cs,g1,an,ath1)
        endif
        if(ffb.gt.(sqrt(ffp**2+0.25)+0.5))then
        an=sqrt(an2)
        call gsy(gamma,ffb,cs,g2,an,ath2)
        endif
      elseif(gamma.ge.gtr)then
        if(ffb.gt.(sqrt(ffp**2+0.25)+0.5))then
        call ssy(gamma,ffb,cs,g1,g2,alpha)
        endif
      endif
c      write(11,*)e1,e2,a1,a2
      e1=e1+g1*em**(-ex)*de*bmag*1.352e-22
      e2=e2+g2*em**(-ex)*de*bmag*1.352e-22
      a1=a1+g1*em**(-ex)*de/ffb/ffb/bmag*1.896e-8*
     &  (ex*gamma*(gamma+1.)+2.*gamma*gamma-1.)/gamma/(gamma**2-1.)
      a2=a2+g2*em**(-ex)*de/ffb/ffb/bmag*1.896e-8*
     &  (ex*gamma*(gamma+1.)+2.*gamma*gamma-1.)/gamma/(gamma**2-1.)
      enddo
c--------end of integration over energy----------------------
      e1d(k)=e1
      e2d(k)=e2
      a1d(k)=a1
      a2d(k)=a2
      if(e2.gt.0.)then
        rccd(k)=(e2-e1)/(e2+e1)
      else
        rccd(k)=0.
      endif
      arg1=anor/2.25e26/th/th*a1
      if(arg1.lt.0.001)then
        phi1(k)=e1*anor/2.25e26
      elseif(arg1.gt.20.)then
        phi1(k)=e1/a1*th*th
      else
        phi1(k)=th*th*e1/a1*(1.-exp(-arg1))
      endif
      arg2=anor/2.25e26/th/th*a2
      if(arg2.lt.0.001)then
        phi2(k)=e2*anor/2.25e26
      elseif(arg2.gt.20.)then
        phi2(k)=e2/a2*th*th
      else
        phi2(k)=th*th*e2/a2*(1.-exp(-arg2))
      endif
c-----calculate Stokes parameters------------------
      phit(k)=phi1(k)+phi2(k)
      q(k)=phi1(k)*(1.-ath1**2)/(1.+ath1**2)+
     &  phi2(k)*(1.-ath2**2)/(1.+ath2**2)
      v(k)=2.*(phi1(k)*ath1/(1.+ath1**2)+phi2(k)*ath2/(1.+ath2**2))
c-----calculate polarization-------------------      
      if(phit(k).gt.0.)then
        rcd(k)=(phi2(k)-phi1(k))/phit(k)
        scd(k)=v(k)/abs(v(k))*sqrt(q(k)**2+v(k)**2)/phit(k)
      else
        rcd(k)=0.
        scd(k)=0.
      endif
c----------------------------------------------
      write(6,*)'ffb= ',ffb
      enddo
c------------end of loop over ffb----------------------------
c
c------------write output------------------------------------
      write(10,100)ex,cs,bmag,anor,theta,j1,j2
      write(10,99)etr,alpha
 100   format(' ex=',f5.2,' cs=',f5.3,' bmag='f5.0,' anor=',e10.2,
     & ' theta=',f5.0,' j1=',i2,' j2=', i3)
99     format(' etr=',f5.2,' alpha= ',f5.3)
      write(10,*)'=================================================='
      write(10,*)'  ffb      emissivities      absorption coefficients
     &    pol'
      write(10,*)'             o mode     x mode     o mode     x mode'
      write(10,*)'           erg/(cm^3 s^1 sr^1 Hz^1)      cm^-1'
      do k=1,kf
      write(10,101)ffbd(k),e1d(k),e2d(k),a1d(k),a2d(k),rccd(k),k
101   format(e11.3,5e11.3,i5)
      enddo
      write(10,*)'=================================================='
      write(10,*)'     f(Hz)        phi_o       phi_x        phi_tot
     &      pol'
      write(10,*)'                        erg cm^-2 s^-1 Hz^-1'    
      do k=1,kf
      ff=ffbd(k)*2.80e6*bmag
      write(10,102)ff,phi1(k)+1.e-23,phi2(k)+1.e-23,
     @   phit(k)+1.e-23,rcd(k),k
102   format(5e13.3,i5)
      enddo
      stop
      end
c--------synchrotron routine---------------------------------
      subroutine ssy(gamma,ffb,cs,g1,g2,alpha)
      dimension xd(34),fd(34)
      data xd/0.001,0.005,0.01,0.025,0.05,0.075,0.1,0.15,0.2,0.25,
     &        0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.,1.2,1.4,1.6,1.8,2.,2.5,
     &        3.,3.5,4.,4.5,5.,6.,7.,8.,9.,10./
      data fd/0.213,0.358,0.445,0.583,0.702,0.722,0.818,0.874,0.904,
     &        0.917,0.919,0.901,0.872,0.832,0.788,0.742,0.694,0.655,
     &        0.566,0.486,0.414,0.354,0.301,0.2,0.13,0.0845,0.0541,
     &        0.0339,0.0214,0.0085,0.0033,0.0013,0.0005,0.00019/
      ss=sqrt(1.-cs*cs)
      ffc=ffb*2./3./ss/gamma/gamma*
     &    (1.+9./4.*(gamma*gamma-1.)/alpha/alpha/ffb/ffb)**1.5
      if(ffc.le.xd(1))then
      f=2.15*ffc**0.333*(1.-0.844*ffc**0.6666)
      elseif(ffc.gt.xd(34))then
      f=1.253*exp(-ffc)*sqrt(ffc)*(1.+55./72/ffc)
      else 
      call sear(xd,fd,ffc,f)
      endif 
      gtot=0.138*ss*f
      g1=gtot/2.
      g2=gtot/2.
      return
      end
c-----------search routine-----------------------------------
      SUBROUTINE SEAR(XC,SGC,X,S)
      DIMENSION XC(34),SGC(34)
      I=1
   12 IF(XC(I).GE.X)GO TO  13
      I=I+1
      GO TO 12
   13 S=(SGC(I)*(X-XC(I-1)  )+SGC(I-1)*(XC(I)   -X ))/(XC(I)-XC(I-1))
      RETURN
      END
c------------gyrosynchrotron routine-------------------------
      subroutine gsy(gamma,ffb,cs,g12,an,ath)
      beta=sqrt(gamma*gamma-1.)/gamma
      ss=sqrt(1.-cs*cs)
      ffc=ffb*2./3./ss/gamma/gamma
      if(ffc.ge.20.)then
      g12=0.
      return
      else
      is1=ffb*gamma*(1.-an*beta*cs)+1
      is2=ffb*gamma*(1.+an*beta*cs)
      sum12=0.
      do i=is1,is2
      cphis=(1.-i/ffb/gamma)/beta/cs/an
      sphis=sqrt(1.-cphis*cphis)
      xs=i*an*beta*ss*sphis/(1.-an*beta*cs*cphis)
      if(ffb.gt.50.)then
      xstr=0.8
      if(gamma.gt.5)xstr=0.9
      if(gamma.gt.10)xstr=0.95
      if(gamma.gt.15)xstr=0.96
      if(xs.lt.xstr*i)goto 10
      endif
      if(ffb.gt.250.)then
      xstr=0.9
      if(gamma.gt.5)xstr=0.92
      if(gamma.gt.10)xstr=0.97
      if(gamma.gt.15)xstr=0.98
      if(xs.lt.xstr*i)goto 10
      endif
      call bes(i,xs,b)
      call bespr(i,xs,bpr)
      f12=(-beta*sphis*bpr+ath*(cs/ss/an-beta*cphis/ss)*b)**2
      s12old=sum12
      sum12=sum12+f12
      if(s12old.gt.0.0)then
      if((sum12-s12old)/s12old.lt.1.e-4)goto 11
      endif
c      write(11,100)gamma,i,xs,f12,sum12
100   format(f10.3,i5,3e12.3)
10    continue
      enddo
11    continue
      g12=sum12/beta/2./cs*ffb/(1.+ath**2)
      return
      endif
      end
c-----index of refraction and polarization coefficient-------
      subroutine refr(ffb,ffp,cs,an1,an2,ath1,ath2)
      ss=sqrt(1.-cs*cs)
      anum=2.*ffp*ffp*(ffp*ffp-ffb*ffb)
      dnum1=+sqrt(ffb**4*ss**4+4.*ffb**2*(ffp**2-ffb**2)**2*cs**2)-
     &     2.*ffb**2*(ffp**2-ffb**2)-ffb**2*ss**2
      dnum2=-sqrt(ffb**4*ss**4+4.*ffb**2*(ffp**2-ffb**2)**2*cs**2)-
     &     2.*ffb**2*(ffp**2-ffb**2)-ffb**2*ss**2
      an1=1.+anum/dnum1
      an2=1.+anum/dnum2
      aknum=2.*ffb*(ffp*ffp-ffb*ffb)*cs
      dknum1=+sqrt(ffb**4*ss**4+4.*ffb**2*(ffp**2-ffb**2)**2*cs**2)-
     &     ffb**2*ss**2
      dknum2=-sqrt(ffb**4*ss**4+4.*ffb**2*(ffp**2-ffb**2)**2*cs**2)-
     &     ffb**2*ss**2
      ath1=-aknum/dknum1
      ath2=-aknum/dknum2
c      write(12,*)ffb,aknum,dknum1,dknum2
      return
      end
c--------------Bessel Function routines----------------------
      subroutine bespr(n,x,bpr)
      n1=n+1
      call bes(n1,x,b1)
      call bes(n,x,b)
      bpr=-b1+n/x*b
      return
      end
      subroutine bes(n,x,b)
      if(n.eq.0)b=bessj0(x)
      if(n.eq.1)b=bessj1(x)
      if(n.ge.2)then
        b=bessj(n,x)
      endif
      return
      end
      FUNCTION BESSJ0(X)
      REAL*8 Y,P1,P2,P3,P4,P5,Q1,Q2,Q3,Q4,Q5,R1,R2,R3,R4,R5,R6,
     *    S1,S2,S3,S4,S5,S6
      DATA P1,P2,P3,P4,P5/1.D0,-.1098628627D-2,.2734510407D-4,
     *    -.2073370639D-5,.2093887211D-6/, Q1,Q2,Q3,Q4,Q5/-.1562499995D-
     *1,
     *    .1430488765D-3,-.6911147651D-5,.7621095161D-6,-.934945152D-7/
      DATA R1,R2,R3,R4,R5,R6/57568490574.D0,-13362590354.D0,651619640.7D
     *0,
     *    -11214424.18D0,77392.33017D0,-184.9052456D0/,
     *    S1,S2,S3,S4,S5,S6/57568490411.D0,1029532985.D0,
     *    9494680.718D0,59272.64853D0,267.8532712D0,1.D0/
      IF(ABS(X).LT.8.)THEN
        Y=X**2
        BESSJ0=(R1+Y*(R2+Y*(R3+Y*(R4+Y*(R5+Y*R6)))))
     *      /(S1+Y*(S2+Y*(S3+Y*(S4+Y*(S5+Y*S6)))))
      ELSE
        AX=ABS(X)
        Z=8./AX
        Y=Z**2
        XX=AX-.785398164
        BESSJ0=SQRT(.636619772/AX)*(COS(XX)*(P1+Y*(P2+Y*(P3+Y*(P4+Y
     *      *P5))))-Z*SIN(XX)*(Q1+Y*(Q2+Y*(Q3+Y*(Q4+Y*Q5)))))
      ENDIF
      RETURN
      END
      FUNCTION BESSJ1(X)
      REAL*8 Y,P1,P2,P3,P4,P5,Q1,Q2,Q3,Q4,Q5,R1,R2,R3,R4,R5,R6,
     *    S1,S2,S3,S4,S5,S6
      DATA R1,R2,R3,R4,R5,R6/72362614232.D0,-7895059235.D0,242396853.1D0
     *,
     *    -2972611.439D0,15704.48260D0,-30.16036606D0/,
     *    S1,S2,S3,S4,S5,S6/144725228442.D0,2300535178.D0,
     *    18583304.74D0,99447.43394D0,376.9991397D0,1.D0/
      DATA P1,P2,P3,P4,P5/1.D0,.183105D-2,-.3516396496D-4,.2457520174D-5
     *,
     *    -.240337019D-6/, Q1,Q2,Q3,Q4,Q5/.04687499995D0,-.2002690873D-3
     *,
     *    .8449199096D-5,-.88228987D-6,.105787412D-6/
      IF(ABS(X).LT.8.)THEN
        Y=X**2
        BESSJ1=X*(R1+Y*(R2+Y*(R3+Y*(R4+Y*(R5+Y*R6)))))
     *      /(S1+Y*(S2+Y*(S3+Y*(S4+Y*(S5+Y*S6)))))
      ELSE
        AX=ABS(X)
        Z=8./AX
        Y=Z**2
        XX=AX-2.356194491
        BESSJ1=SQRT(.636619772/AX)*(COS(XX)*(P1+Y*(P2+Y*(P3+Y*(P4+Y
     *      *P5))))-Z*SIN(XX)*(Q1+Y*(Q2+Y*(Q3+Y*(Q4+Y*Q5)))))
     *      *SIGN(1.,X)
      ENDIF
      RETURN
      END
      FUNCTION BESSJ(N,X)
      PARAMETER (IACC=40,BIGNO=1.E10,BIGNI=1.E-10)
      IF(N.LT.2)PAUSE 'bad argument N in BESSJ'
      TOX=2./X
      IF(X.GT.FLOAT(N))THEN
        BJM=BESSJ0(X)
        BJ=BESSJ1(X)
        DO 11 J=1,N-1
          BJP=J*TOX*BJ-BJM
          BJM=BJ
          BJ=BJP
11      CONTINUE
        BESSJ=BJ
      ELSE
        M=2*((N+INT(SQRT(FLOAT(IACC*N))))/2)
        BESSJ=0.
        JSUM=0
        SUM=0.
        BJP=0.
        BJ=1.
        DO 12 J=M,1,-1
          BJM=J*TOX*BJ-BJP
          BJP=BJ
          BJ=BJM
          IF(ABS(BJ).GT.BIGNO)THEN
            BJ=BJ*BIGNI
            BJP=BJP*BIGNI
            BESSJ=BESSJ*BIGNI
            SUM=SUM*BIGNI
          ENDIF
          IF(JSUM.NE.0)SUM=SUM+BJ
          JSUM=1-JSUM
          IF(J.EQ.N)BESSJ=BJP
12      CONTINUE
        SUM=2.*SUM-BJ
        BESSJ=BESSJ/SUM
      ENDIF
      RETURN
      END