      SUBROUTINE BSSLI (MO, A, IN, W)
C     ******************************************************************
C     FORTRAN SUBROUTINE FOR MODIFIED BESSEL FUNCTION OF INTEGRAL ORDER
C     ******************************************************************
C     MO = MODE OF OPERATION
C     A  = ARGUMENT (COMPLEX NUMBER)
C     IN = ORDER (INTEGER)
C     W  = FUNCTION OF FIRST KIND (COMPLEX NUMBER)
C     -------------------
      COMPLEX A, W
      DIMENSION AZ(2), FI(2)
      DIMENSION CD(30), CE(30)
      DIMENSION QZ(2), RZ(2), SZ(2), ZR(2)
      DIMENSION TS(2), TM(2), RM(4), SM(4), AQ(2), QF(2)
      DATA CD(1) / 0.00000000000000E00/,  CD(2) /-1.64899505142212E-2/,
     1     CD(3) /-7.18621880068536E-2/,  CD(4) /-1.67086878124866E-1/,
     2     CD(5) /-3.02582250219469E-1/,  CD(6) /-4.80613945245927E-1/,
     3     CD(7) /-7.07075239357898E-1/,  CD(8) /-9.92995790539516E-1/,
     4     CD(9) /-1.35583925612592E00/,  CD(10)/-1.82105907899132E00/,
     5     CD(11)/-2.42482175310879E00/,  CD(12)/-3.21956655708750E00/,
     6     CD(13)/-4.28658077248384E00/,  CD(14)/-5.77022816798128E00/,
     7     CD(15)/-8.01371260952526E00/
      DATA CD(16)/ 0.00000000000000E00/,  CD(17)/-5.57742429879505E-3/,
     1     CD(18)/-4.99112944172476E-2/,  CD(19)/-1.37440911652397E-1/,
     2     CD(20)/-2.67233784710566E-1/,  CD(21)/-4.40380166808682E-1/,
     3     CD(22)/-6.61813614872541E-1/,  CD(23)/-9.41861077665017E-1/,
     4     CD(24)/-1.29754130468326E00/,  CD(25)/-1.75407696719816E00/,
     5     CD(26)/-2.34755299882276E00/,  CD(27)/-3.13041332689196E00/,
     6     CD(28)/-4.18397120563729E00/,  CD(29)/-5.65251799214994E00/,
     7     CD(30)/-7.87863959810677E00/
      DATA CE(1) / 0.00000000000000E00/,  CE(2) /-4.80942336387447E-3/,
     1     CE(3) /-1.31366200347759E-2/,  CE(4) /-1.94843834008458E-2/,
     2     CE(5) /-2.19948900032003E-2/,  CE(6) /-2.09396625676519E-2/,
     3     CE(7) /-1.74600268458650E-2/,  CE(8) /-1.27937813362085E-2/,
     4     CE(9) /-8.05234421796592E-3/,  CE(10)/-4.15817375002760E-3/,
     5     CE(11)/-1.64317738747922E-3/,  CE(12)/-4.49175585314709E-4/,
     6     CE(13)/-7.28594765574007E-5/,  CE(14)/-5.38265230658285E-6/,
     7     CE(15)/-9.93779048036289E-8/
      DATA CE(16)/ 0.00000000000000E00/,  CE(17)/ 7.53805779200591E-2/,
     1     CE(18)/ 7.12293537403464E-2/,  CE(19)/ 6.33116224228200E-2/,
     2     CE(20)/ 5.28240264523301E-2/,  CE(21)/ 4.13305359441492E-2/,
     3     CE(22)/ 3.01350573947510E-2/,  CE(23)/ 2.01043439592720E-2/,
     4     CE(24)/ 1.18552223068074E-2/,  CE(25)/ 5.86055510956010E-3/,
     5     CE(26)/ 2.25465148267325E-3/,  CE(27)/ 6.08173041536336E-4/,
     6     CE(28)/ 9.84215550625747E-5/,  CE(29)/ 7.32139093038089E-6/,
     7     CE(30)/ 1.37279667384666E-7/
C     -------------------
      AZ(1)=REAL(A)
      AZ(2)=AIMAG(A)
      ZS=AZ(1)*AZ(1)+AZ(2)*AZ(2)
      ZM=SQRT(ZS)
      PN=IABS(IN)
      SN=+1.0
      IF(AZ(1))002,003,003
  002 QZ(1)=-AZ(1)
      QZ(2)=-AZ(2)
      IF(IN.EQ.IN/2*2)GO TO 004
      SN=-1.0
      GO TO 004
  003 QZ(1)=AZ(1)
      QZ(2)=AZ(2)
  004 IF(ZM.LE.17.5+0.5*PN*PN)GO TO 005
      QN=PN
      GO TO 011
  005 QN=0.5*ZM-0.5*ABS(QZ(1))+0.5*ABS(0.5*ZM-ABS(QZ(1)))
      IF(PN.LE.QN)GO TO 006
      QN=+AINT(0.0625*ZS)
      IF(PN.LE.QN)GO TO 039
      QN=PN
      GO TO 039
  006 IF(ZM.LE.17.5)GO TO 007
      QN=+AINT(SQRT(2.0*(ZM-17.5)))
      GO TO 011
  007 IF(ZS-1.0)009,008,008
  008 IF(-ABS(AZ(1))+0.096*AZ(2)*AZ(2))009,010,010
  009 QN=AINT(0.0625*ZS)
      IF(PN.LE.QN)GO TO 039
      QN=PN
      GO TO 039
  010 QN=0.0
  011 SZ(1)=QZ(1)
      SZ(2)=QZ(2)
      QM=SN*0.398942280401433
      ZR(1)=SQRT(SZ(1)+ZM)
      ZR(2)=SZ(2)/ZR(1)
      ZR(1)=0.707106781186548*ZR(1)
      ZR(2)=0.707106781186548*ZR(2)
      QF(1)=+QM*ZR(1)/ZM
      QF(2)=-QM*ZR(2)/ZM
      IF(ZM.LE.17.5)GO TO 017
  012 RZ(1)=+0.5*QZ(1)/ZS
      RZ(2)=-0.5*QZ(2)/ZS
      AN=QN*QN-0.25
      SM(1)=0.0
      SM(2)=0.0
      SM(3)=0.0
      SM(4)=0.0
      TM(1)=1.0
      TM(2)=0.0
      PM=0.0
      GO TO 014
  013 AN=AN-2.0*PM
      PM=PM+1.0
      TS(1)=TM(1)*RZ(1)-TM(2)*RZ(2)
      TS(2)=TM(1)*RZ(2)+TM(2)*RZ(1)
      TM(1)=AN*TS(1)/PM
      TM(2)=AN*TS(2)/PM
  014 SM(1)=SM(1)+TM(1)
      SM(2)=SM(2)+TM(2)
      AN=AN-2.0*PM
      PM=PM+1.0
      TS(1)=TM(1)*RZ(1)-TM(2)*RZ(2)
      TS(2)=TM(1)*RZ(2)+TM(2)*RZ(1)
      TM(1)=AN*TS(1)/PM
      TM(2)=AN*TS(2)/PM
      IF(ABS(SM(3))+ABS(TM(1)).NE.ABS(SM(3)))GO TO 015
      IF(ABS(SM(4))+ABS(TM(2)).EQ.ABS(SM(4)))GO TO 016
  015 SM(3)=SM(3)+TM(1)
      SM(4)=SM(4)+TM(2)
      IF(PM.LT.35.0)GO TO 013
  016 TS(1)=SM(1)+SM(3)
      TS(2)=SM(2)+SM(4)
      SM(1)=SM(1)-SM(3)
      SM(2)=SM(2)-SM(4)
      SM(3)=TS(1)
      SM(4)=TS(2)
      GO TO 019
  017 SM(1)=1.0
      SM(2)=0.0
      SM(3)=1.0
      SM(4)=0.0
      M=15.0*QN+2.0
      N=15.0*QN+15.0
      DO 018 I=M,N
      TS(1)=-QZ(1)-CD(I)
      TS(2)=-QZ(2)
      SS=TS(1)*TS(1)+TS(2)*TS(2)
      TM(1)=+CE(I)*TS(1)/SS
      TM(2)=-CE(I)*TS(2)/SS
      SM(1)=SM(1)+TM(1)
      SM(2)=SM(2)+TM(2)
      TS(1)=QZ(1)-CD(I)
      TS(2)=QZ(2)
      SS=TS(1)*TS(1)+TS(2)*TS(2)
      TM(1)=+CE(I)*TS(1)/SS
      TM(2)=-CE(I)*TS(2)/SS
      SM(3)=SM(3)+TM(1)
      SM(4)=SM(4)+TM(2)
  018 CONTINUE
  019 RM(1)=SM(1)
      RM(2)=SM(2)
      IF(QZ(1).GE.17.5)GO TO 023
      AQ(1)=-2.0*QZ(1)
      IF(QZ(2))020,021,021
  020 AQ(2)=-2.0*QZ(2)-3.14159265358979*(QN+0.5)
      GO TO 022
  021 AQ(2)=-2.0*QZ(2)+3.14159265358979*(QN+0.5)
  022 QM=EXP(AQ(1))
      TS(1)=QM*COS(AQ(2))
      TS(2)=QM*SIN(AQ(2))
      RM(1)=RM(1)+TS(1)*SM(3)-TS(2)*SM(4)
      RM(2)=RM(2)+TS(1)*SM(4)+TS(2)*SM(3)
  023 IF(QN.EQ.PN)GO TO 037
      RM(3)=RM(1)
      RM(4)=RM(2)
      QN=QN+1.0
      IF(ZM.LE.17.5)GO TO 029
  024 AN=QN*QN-0.25
      SM(1)=0.0
      SM(2)=0.0
      SM(3)=0.0
      SM(4)=0.0
      TM(1)=1.0
      TM(2)=0.0
      PM=0.0
      GO TO 026
  025 AN=AN-2.0*PM
      PM=PM+1.0
      TS(1)=TM(1)*RZ(1)-TM(2)*RZ(2)
      TS(2)=TM(1)*RZ(2)+TM(2)*RZ(1)
      TM(1)=AN*TS(1)/PM
      TM(2)=AN*TS(2)/PM
  026 SM(1)=SM(1)+TM(1)
      SM(2)=SM(2)+TM(2)
      AN=AN-2.0*PM
      PM=PM+1.0
      TS(1)=TM(1)*RZ(1)-TM(2)*RZ(2)
      TS(2)=TM(1)*RZ(2)+TM(2)*RZ(1)
      TM(1)=AN*TS(1)/PM
      TM(2)=AN*TS(2)/PM
      IF(ABS(SM(3))+ABS(TM(1)).NE.ABS(SM(3)))GO TO 027
      IF(ABS(SM(4))+ABS(TM(2)).EQ.ABS(SM(4)))GO TO 028
  027 SM(3)=SM(3)+TM(1)
      SM(4)=SM(4)+TM(2)
      IF(PM.LT.35.0)GO TO 025
  028 TS(1)=SM(1)+SM(3)
      TS(2)=SM(2)+SM(4)
      SM(1)=SM(1)-SM(3)
      SM(2)=SM(2)-SM(4)
      SM(3)=TS(1)
      SM(4)=TS(2)
      GO TO 031
  029 SM(1)=1.0
      SM(2)=0.0
      SM(3)=1.0
      SM(4)=0.0
      M=15.0*QN+2.0
      N=15.0*QN+15.0
      DO 030 I=M,N
      TS(1)=-QZ(1)-CD(I)
      TS(2)=-QZ(2)
      SS=TS(1)*TS(1)+TS(2)*TS(2)
      TM(1)=+CE(I)*TS(1)/SS
      TM(2)=-CE(I)*TS(2)/SS
      SM(1)=SM(1)+TM(1)
      SM(2)=SM(2)+TM(2)
      TS(1)=+QZ(1)-CD(I)
      TS(2)=+QZ(2)
      SS=TS(1)*TS(1)+TS(2)*TS(2)
      TM(1)=+CE(I)*TS(1)/SS
      TM(2)=-CE(I)*TS(2)/SS
      SM(3)=SM(3)+TM(1)
      SM(4)=SM(4)+TM(2)
  030 CONTINUE
  031 RM(1)=SM(1)
      RM(2)=SM(2)
      IF(QZ(1).GE.17.5)GO TO 036
      AQ(1)=-2.0*QZ(1)
      IF(QZ(2))032,033,033
  032 AQ(2)=-2.0*QZ(2)-3.14159265358979*(QN+0.5)
      GO TO 034
  033 AQ(2)=-2.0*QZ(2)+3.14159265358979*(QN+0.5)
  034 QM=EXP(AQ(1))
      TS(1)=QM*COS(AQ(2))
      TS(2)=QM*SIN(AQ(2))
      RM(1)=RM(1)+TS(1)*SM(3)-TS(2)*SM(4)
      RM(2)=RM(2)+TS(1)*SM(4)+TS(2)*SM(3)
      GO TO 036
  035 TM(1)=-2.0*QN*QZ(1)/ZS
      TM(2)=+2.0*QN*QZ(2)/ZS
      TS(1)=TM(1)*RM(1)-TM(2)*RM(2)+RM(3)
      TS(2)=TM(1)*RM(2)+TM(2)*RM(1)+RM(4)
      RM(3)=RM(1)
      RM(4)=RM(2)
      RM(1)=TS(1)
      RM(2)=TS(2)
      QN=QN+1.0
  036 IF(QN.LT.PN)GO TO 035
  037 IF(MO.NE.0)GO TO 038
      QM=EXP(QZ(1))
      TM(1)=QM*COS(QZ(2))
      TM(2)=QM*SIN(QZ(2))
      TS(1)=TM(1)*RM(1)-TM(2)*RM(2)
      TS(2)=TM(1)*RM(2)+TM(2)*RM(1)
      RM(1)=TS(1)
      RM(2)=TS(2)
  038 FI(1)=QF(1)*RM(1)-QF(2)*RM(2)
      FI(2)=QF(1)*RM(2)+QF(2)*RM(1)
      W=CMPLX(FI(1),FI(2))
      RETURN
  039 SZ(1)=0.25*(QZ(1)*QZ(1)-QZ(2)*QZ(2))
      SZ(2)=0.5*QZ(1)*QZ(2)
      AN=QN
      SM(1)=0.0
      SM(2)=0.0
      SM(3)=0.0
      SM(4)=0.0
      TM(1)=1.0
      TM(2)=0.0
      PM=0.0
  040 AN=AN+1.0
      TS(1)=TM(1)/AN
      TS(2)=TM(2)/AN
      SM(3)=SM(3)+TS(1)
      SM(4)=SM(4)+TS(2)
      TM(1)=TS(1)*SZ(1)-TS(2)*SZ(2)
      TM(2)=TS(1)*SZ(2)+TS(2)*SZ(1)
      PM=PM+1.0
      TM(1)=TM(1)/PM
      TM(2)=TM(2)/PM
      IF(ABS(SM(1))+ABS(TM(1)).NE.ABS(SM(1)))GO TO 041
      IF(ABS(SM(2))+ABS(TM(2)).EQ.ABS(SM(2)))GO TO 042
  041 SM(1)=SM(1)+TM(1)
      SM(2)=SM(2)+TM(2)
      GO TO 040
  042 SM(1)=SM(1)+1.0
      AN=QN+1.0
      SM(3)=AN*SM(3)
      SM(4)=AN*SM(4)
      GO TO 044
  043 AN=QN*(QN+1.0)
      TM(1)=SZ(1)/AN
      TM(2)=SZ(2)/AN
      TS(1)=+TM(1)*SM(3)-TM(2)*SM(4)
      TS(2)=+TM(1)*SM(4)+TM(2)*SM(3)
      SM(3)=SM(1)
      SM(4)=SM(2)
      SM(1)=SM(1)+TS(1)
      SM(2)=SM(2)+TS(2)
      QN=QN-1.0
  044 IF(QN.GT.PN)GO TO 043
      QF(1)=SN
      QF(2)=0.0
      QN=0.0
      GO TO 046
  045 QN=QN+1.0
      TM(1)=QF(1)*QZ(1)-QF(2)*QZ(2)
      TM(2)=QF(1)*QZ(2)+QF(2)*QZ(1)
      QF(1)=0.5*TM(1)/QN
      QF(2)=0.5*TM(2)/QN
  046 IF(QN.LT.PN)GO TO 045
      IF(MO.EQ.0)GO TO 047
      QM=EXP(-QZ(1))
      TM(1)=QM*COS(-QZ(2))
      TM(2)=QM*SIN(-QZ(2))
      TS(1)=TM(1)*QF(1)-TM(2)*QF(2)
      TS(2)=TM(1)*QF(2)+TM(2)*QF(1)
      QF(1)=TS(1)
      QF(2)=TS(2)
  047 FI(1)=QF(1)*SM(1)-QF(2)*SM(2)
      FI(2)=QF(1)*SM(2)+QF(2)*SM(1)
      W=CMPLX(FI(1),FI(2))
      RETURN
      END
