      Subroutine light(phs,xincl,xh,xc,yh,yc,n1,n2,sumhot,sumkul,rv,grx,   15650
     $gry,grz,rvq,grxq,gryq,grzq,mmsave,theta,rho,aa,bb,slump1,slump2,     15660
     $somhot,somkul,d,wl,snth,csth,snfi,csfi,tld,gmag1,gmag2)              15670
c  Version of July 27, 1992                          
      DIMENSION RV(*),GRX(*),GRY(*),GRZ(*),RVQ(*),GRXQ(*),GRYQ(*),GRZQ(*   15690
     $),SLUMP1(*),SLUMP2(*),MMSAVE(*),THETA(*),RHO(*),AA(*),BB(*)          15700
      DIMENSION SNTH(*),CSTH(*),SNFI(*),CSFI(*),tld(*),gmag1(*),gmag2(*)   15710
      COMMON X1                                                            15720
      COMMON /KFAC/ KFF1,KFF2                                              15730
      COMMON /NSPT/ NSP1,NSP2,IFAT1,IFAT2                                  15740
      common /invar/ khdum,ipbdum,irtedm,nrefdm,irv1dm,irv2dm,mrefdm       15750
     $,ifs1dm,ifs2dm,icr1dm,icr2dm,ld                                      15760
      common /qrm/ rm
c     Added by C. Hummel
      integer*4 maxlc
      parameter (maxlc=1000)
      real*8 xlc(maxlc),ylc(maxlc),flc(maxlc)
      integer*4 clc(maxlc),ilc
      common /imlcdata1/clc,ilc
      common /imlcdata2/xlc,ylc,flc
      ilc=0

      rmf=d*rm/(1.+rm)
      PHA=PHS*6.283185                                                     15770
      K=6                                                                  15780
      KK=K+1                                                               15790
      XLUMP=14384./WL                                                      15800
      XINC=XINCL*.0174533                                                  15810
      L=1                                                                  15820
      TEST=(PHS-.5)**2                                                     15830
      TESTS=(TEST-.071525)**2                                              15840
      SINI=SIN(XINC)                                                       15850
      COSPH=COS(PHA)                                                       15860
      SINPH=SIN(PHA)                                                       15870
      SINSQ=SINPH**2                                                       15880
      COSI=COS(XINC)                                                       15890
      NP1=N1+1                                                             15900
      NP2=N1+N2+2                                                          15910
      LLL1=MMSAVE(NP1)                                                     15920
      LLL2=MMSAVE(NP2)                                                     15930
      NPP2=NP2-1                                                           15940
      LL1=MMSAVE(N1)+1                                                     15950
      LL2=MMSAVE(NPP2)+1                                                   15960
      LLLL1=(LL1+LLL1)/2                                                   15970
      LLLL2=(LL2+LLL2)/2                                                   15980
      SINSQE=0.                                                            15990
      IF(SINI.GT.0.) SINSQE=((1.02*(RV(LLL1)+RVQ(LLL2))/D)**2-COSI**2)/    16000
     $SINI**2                                                              16010
      CICP=COSI*COSPH                                                      16020
      CISP=COSI*SINPH                                                      16030
      XLOS=COSPH*SINI                                                      16040
      YLOS=-SINPH*SINI                                                     16050
      ZLOS=COSI                                                            16060
      SUM=0.                                                               16070
      SOM=0.                                                               16080
      IF(TEST.LE..0625)GOTO 18                                             16090
      COMP=-1.                                                             16100
      CMP=1.                                                               16110
      COMPP=1.                                                             16120
      KOMP=2                                                               16130
      NSPOT=NSP2                                                           16140
      IFAT=IFAT2                                                           16150
      CMPP=0.                                                              16160
      X=XC                                                                 16170
      y=yc                                                                 16180
      EN=N2                                                                16190
      NPH=N2                                                               16200
      NP=2*N2                                                              16210
      GOTO 28                                                              16220
   18 X=XH                                                                 16230
      y=yh                                                                 16240
      COMP=1.                                                              16250
      KOMP=1                                                               16260
      NSPOT=NSP1                                                           16270
      IFAT=IFAT1                                                           16280
      CMP=0.                                                               16290
      COMPP=-1.                                                            16300
      CMPP=1.                                                              16310
      EN=N1                                                                16320
      NPH=N1                                                               16330
      NP=2*N1                                                              16340
   28 DELTH=1.570796/EN                                                    16350
      AR=CMPP*RV(LLLL1)+CMP*RVQ(LLLL2)                                     16360
      BR=CMPP*RV(1)+CMP*RVQ(1)                                             16370
      ASQ=AR*AR                                                            16380
      BSQ=BR*BR                                                            16390
      AB=AR*BR                                                             16400
      ASBS=ASQ-BSQ                                                         16410
      KF=(2-KOMP)*KFF1+(KOMP-1)*KFF2                                       16420
      CMPPD=CMPP*D                                                         16430
      CMPD=CMP*D                                                           16440
      NPP=NP+1                                                             16450
      TEMF=1.                                                              16460
      DO 36 I=1,NP                                                         16470
      IF(I.GT.NPH)GOTO 54                                                  16480
      UPDOWN=1.                                                            16490
      IK=I                                                                 16500
      GOTO 55                                                              16510
   54 UPDOWN=-1.                                                           16520
      IK=NPP-I                                                             16530
   55 CONTINUE                                                             16540
      IPN1=IK+(KOMP-1)*N1                                                  16550
      SINTH=SNTH(IPN1)                                                     16560
      COSTH=CSTH(IPN1)*UPDOWN                                              16570
      EM=SINTH*EN*1.3                                                      16580
      MM=EM+1.                                                             16590
      XM=MM                                                                16600
      MH=MM                                                                16610
      MM=2*MM                                                              16620
      DELFI=3.141593/XM                                                    16630
      IP=(KOMP-1)*NP1+IK                                                   16640
      IY=MMSAVE(IP)+1                                                      16650
      IF(TEST.LE..0625)GOTO 19                                             16660
      GX=GRXQ(IY)                                                          16670
      GY=-GRYQ(IY)                                                         16680
      GZ=UPDOWN*GRZQ(IY)                                                   16690
      grmag=gmag2(iy)                                                      16700
      GOTO 29                                                              16710
   19 GX=GRX(IY)                                                           16720
      GY=-GRY(IY)                                                          16730
      GZ=UPDOWN*GRZ(IY)                                                    16740
      grmag=gmag1(iy)                                                      16750
   29 COSSAV=(XLOS*GX+YLOS*GY+ZLOS*GZ)/GRMAG                               16760
      SUMJ=0.                                                              16770
      SOMJ=0.                                                              16780
      MPP=MM+1                                                             16790
      IY=IY-1                                                              16800
      DO 26 J=1,MM                                                         16810
      IF(J.GT.MH) GOTO 58                                                  16820
      RTLEFT=1.                                                            16830
      JK=J                                                                 16840
      GOTO 59                                                              16850
   58 RTLEFT=-1.                                                           16860
      JK=MPP-J                                                             16870
   59 CONTINUE                                                             16880
      IX=IY+JK                                                             16890
      IS=IX+(KOMP-1)*LLL1                                                  16900
      SINFI=SNFI(IS)*RTLEFT                                                16910
      COSFI=CSFI(IS)                                                       16920
      STSF=SINTH*SINFI                                                     16930
      STCF=SINTH*COSFI                                                     16940
      IF(TEST.LE..0625)GOTO 39                                             16950
      IF(RVQ(IX).EQ.-1.) GOTO 26                                           16960
      GX=GRXQ(IX)                                                          16970
      GY=RTLEFT*GRYQ(IX)                                                   16980
      GZ=UPDOWN*GRZQ(IX)                                                   16990
      R=RVQ(IX)                                                            17000
      grmag=gmag2(ix)                                                      17010
      GOTO 49                                                              17020
   39 IF(RV(IX).EQ.-1.) GOTO 26                                            17030
      GX=GRX(IX)                                                           17040
      GY=RTLEFT*GRY(IX)                                                    17050
      GZ=UPDOWN*GRZ(IX)                                                    17060
      R=RV(IX)                                                             17070
      grmag=gmag1(ix)                                                      17080
   49 COSGAM=(XLOS*GX+YLOS*GY+ZLOS*GZ)/GRMAG                               17090
      IF(TESTS.LT.2.2562E-3) GOTO 170                                      17100
      IF((STCF*R).GT.(.7*(CMP+COMP*X1))) GOTO 129                          17110
  170 PROD=COSSAV*COSGAM                                                   17120
      IF(PROD.GT.0.) GOTO 22                                               17130
      COSSAV=-COSSAV                                                       17140
      ZZ=R*COSTH                                                           17150
      YY=R*COMP*STSF                                                       17160
      XX=CMPD+COMP*STCF*R                                                  17170
      YSKY=XX*SINPH+YY*COSPH-cmpd*SINPH                                    17180
      ZSKY=-XX*CICP+yy*CISP+ZZ*SINI+CMPD*CICP                              17190
      RHO(L)=SQRT(YSKY**2+ZSKY**2)                                         17200
      THETA(L)=ARSIN(ZSKY/RHO(L))                                          17210
      IF(YSKY.LT.0.) GOTO 92                                               17220
      THETA(L)=6.283185+THETA(L)                                           17230
      GOTO 93                                                              17240
   92 THETA(L)=3.141593-THETA(L)                                           17250
   93 IF (THETA(L).GE.6.283185) THETA(L)=THETA(L)-6.283185                 17260
      L=L+1                                                                17270
      GOTO 27                                                              17280
   22 COSSAV=COSGAM                                                        17290
      GOTO 27                                                              17300
  129 COSSAV=COSGAM                                                        17310
      IF(KF.LE.0) GOTO 27                                                  17320
      ZZ=R*COSTH                                                           17330
      YY=R*COMP*STSF                                                       17340
      XX=CMPD+COMP*STCF*R                                                  17350
      YSKY=XX*SINPH+YY*COSPH-cmpd*SINPH                                    17360
      ZSKY=-XX*CICP+YY*CISP+ZZ*SINI+CMPD*CICP                              17370
      RPT=SQRT(YSKY**2+ZSKY**2)                                            17380
      RADTST=AB/SQRT(BSQ+ASBS*(ZSKY/RPT)**2)                               17390
      IF(RPT.LE.RADTST) GOTO 26                                            17400
   27 IF(COSGAM.GE.0.) GOTO 26                                             17410
      COSGAM=-COSGAM                                                       17420
      DARKEN=1.-X+X*COSGAM                                                 17430
      if(ld.eq.1) goto 141                                                 17440
      if(cosgam.eq.0.) goto 141                                            17450
      darken=darken-y*cosgam*alog(cosgam)                                  17460
      if(darken.lt.0.) darken=0.                                           17470
  141 continue                                                             17480
      ATRATN=1.                                                            17490
      ATRATO=1.                                                            17500
      CORFAC=1.                                                            17510
      IF(NSPOT.EQ.0) GOTO 640                                              17520
      CALL SPOT(KOMP,NSPOT,SINTH,COSTH,SINFI,COSFI,TEMF)                   17530
      IF(TEMF.EQ.1.) GOTO 640                                              17540
      TSP=TLD(IS)*TEMF                                                     17550
      IF(IFAT.EQ.0) GOTO 641                                               17560
      CALL ATM(TLD(IS),WL,ATRATO)                                          17570
      CALL ATM(TSP,WL,ATRATN)                                              17580
  641 CORFAC=ATRATN*(EXP(XLUMP/TLD(IS))-1.)/(ATRATO*(EXP(XLUMP/TSP)-1.))   17590
  640 CONTINUE                                                             17600
      DIF=COSGAM*DARKEN*CORFAC*(CMP*SLUMP2(IX)+CMPP*SLUMP1(IX))            17610
      dift=dif*delfi*delth
      xx=cmpd+comp*stcf*r
      yy=r*comp*stsf
      zz=r*costh
      xcm=xx-rmf
      ys=xcm*sinph+yy*cosph
      zs=-xcm*cicp+yy*cisp+zz*sini

ch    write(6,888) komp,phs,ys,zs,dift
      ilc=ilc+1         ! added by C. Hummel
      xlc(ilc)=dble(ys)
      ylc(ilc)=dble(zs)
      flc(ilc)=dble(dift)
      clc(ilc)=komp

  888 format(1x,i1,4f16.7)
      v=-r*(STCF*YLOS-stsf*XLOS)*COMP                                      17620
      DIFF=DIF*V                                                           17630
      SOMJ=SOMJ+DIFF                                                       17640
   45 SUMJ=SUMJ+DIF                                                        17650
   26 CONTINUE                                                             17660
      SOMJ=SOMJ*DELFI                                                      17670
      SUMJ=SUMJ*DELFI                                                      17680
      SOM=SOM+SOMJ                                                         17690
   36 SUM=SUM+SUMJ                                                         17700
      IF(SINSQ.GE.SINSQE) GOTO 75                                          17710
      L=L-1                                                                17720
      CALL SORTER (L,THETA,RHO)                                            17730
      CALL FOUR (L,K,THETA,RHO,AA,BB)                                      17740
   75 IF(TEST.LE..0625)GOTO 118                                            17750
      SUMKUL=SUM*DELTH                                                     17760
      SOMKUL=SOM*DELTH                                                     17770
      X=XH                                                                 17780
      y=yh                                                                 17790
      KOMP=1                                                               17800
      NSPOT=NSP1                                                           17810
      IFAT=IFAT1                                                           17820
      EN=N1                                                                17830
      SAFTY=2.6*RV(LLL1)/EN                                                17840
      RMAX=RVQ(LLL2)+SAFTY                                                 17850
      RMIN=RVQ(1)-SAFTY                                                    17860
      NPH=N1                                                               17870
      NP=2*N1                                                              17880
      GOTO 128                                                             17890
  118 X=XC                                                                 17900
      y=yc                                                                 17910
      KOMP=2                                                               17920
      NSPOT=NSP2                                                           17930
      IFAT=IFAT2                                                           17940
      SUMHOT=SUM*DELTH                                                     17950
      SOMHOT=SOM*DELTH                                                     17960
      EN=N2                                                                17970
      SAFTY=2.6*RVQ(LLL2)/EN                                               17980
      RMAX=RV(LLL1)+SAFTY                                                  17990
      RMIN=RV(1)-SAFTY                                                     18000
      NPH=N2                                                               18010
      NP=2*N2                                                              18020
  128 DELTH=1.570796/EN                                                    18030
      SOM=0.                                                               18040
      SUM=0.                                                               18050
      NPP=NP+1                                                             18060
      TEMF=1.                                                              18070
      DO 136 I=1,NP                                                        18080
      IF(I.GT.NPH) GOTO 154                                                18090
      UPDOWN=1.                                                            18100
      IK=I                                                                 18110
      GOTO 155                                                             18120
  154 UPDOWN=-1.                                                           18130
      IK=NPP-I                                                             18140
  155 CONTINUE                                                             18150
      IPN1=IK+(KOMP-1)*N1                                                  18160
      SINTH=SNTH(IPN1)                                                     18170
      COSTH=CSTH(IPN1)*UPDOWN                                              18180
      EM=SINTH*EN*1.3                                                      18190
      MM=EM+1.                                                             18200
      XM=MM                                                                18210
      MH=MM                                                                18220
      MM=2*MM                                                              18230
      DELFI=3.141593/XM                                                    18240
      SOMJ=0.                                                              18250
      SUMJ=0.                                                              18260
      SIGN=0.                                                              18270
      DRHO=1.                                                              18280
      MPP=MM+1                                                             18290
      DO 126 J=1,MM                                                        18300
      IF(J.GT.MH) GOTO 158                                                 18310
      RTLEFT=1.                                                            18320
      JK=J                                                                 18330
      GOTO 159                                                             18340
  158 RTLEFT=-1.                                                           18350
      JK=MPP-J                                                             18360
  159 CONTINUE                                                             18370
      IP=(KOMP-1)*NP1+IK                                                   18380
      IX=MMSAVE(IP)+JK                                                     18390
      IS=IX+LLL1*(KOMP-1)                                                  18400
      SINFI=SNFI(IS)*RTLEFT                                                18410
      COSFI=CSFI(IS)                                                       18420
      STSF=SINTH*SINFI                                                     18430
      STCF=SINTH*COSFI                                                     18440
      IF(TEST.LE..0625)GOTO 139                                            18450
      IF(RV(IX).EQ.-1.) GOTO 126                                           18460
      GX=GRX(IX)                                                           18470
      GY=RTLEFT*GRY(IX)                                                    18480
      GZ=UPDOWN*GRZ(IX)                                                    18490
      R=RV(IX)                                                             18500
      grmag=gmag1(ix)                                                      18510
      GOTO 149                                                             18520
  139 IF(RVQ(IX).EQ.-1.) GOTO 126                                          18530
      GX=GRXQ(IX)                                                          18540
      GY=RTLEFT*GRYQ(IX)                                                   18550
      GZ=UPDOWN*GRZQ(IX)                                                   18560
      R=RVQ(IX)                                                            18570
      grmag=gmag2(ix)                                                      18580
  149 COSGAM=(XLOS*GX+YLOS*GY+ZLOS*GZ)/GRMAG                               18590
      IF(COSGAM.LT.0.) GOTO 104                                            18600
      SIGN=0.                                                              18610
      OLSIGN=0.                                                            18620
      GOTO 126                                                             18630
  104 COSGAM=-COSGAM                                                       18640
      ZZ=R*COSTH                                                           18650
      YY=R*COMPP*STSF                                                      18660
      XX=CMPPD+COMPP*STCF*R                                                18670
      DARKEN=1.-X+X*COSGAM                                                 18680
      if(ld.eq.1) goto 142                                                 18690
      if(cosgam.eq.0.) goto 142                                            18700
      darken=darken-y*cosgam*alog(cosgam)                                  18710
      if(darken.lt.0.) darken=0.                                           18720
  142 continue                                                             18730
      OLDIF=DIF                                                            18740
      ATRATN=1.                                                            18750
      ATRATO=1.                                                            18760
      CORFAC=1.                                                            18770
      IF(NSPOT.EQ.0) GOTO 660                                              18780
      CALL SPOT(KOMP,NSPOT,SINTH,COSTH,SINFI,COSFI,TEMF)                   18790
      IF(TEMF.EQ.1.) GOTO 660                                              18800
      TSP=TLD(IS)*TEMF                                                     18810
      IF(IFAT.EQ.0) GOTO 661                                               18820
      CALL ATM(TLD(IS),WL,ATRATO)                                          18830
      CALL ATM(TSP,WL,ATRATN)                                              18840
  661 CORFAC=ATRATN*(EXP(XLUMP/TLD(IS))-1.)/(ATRATO*(EXP(XLUMP/TSP)-1.))   18850
  660 CONTINUE                                                             18860
      DIF=COSGAM*DARKEN*CORFAC*(CMPP*SLUMP2(IX)+CMP*SLUMP1(IX))            18870
      v=R*(STCF*YLOS-STSF*XLOS)*COMP                                       18880
      DIFF=DIF*V                                                           18890
      IF(SINSQ.GT.SINSQE) GOTO 63                                          18900
      OLSIGN=SIGN                                                          18910
      OLDRHO=DRHO                                                          18920
      YSKY=XX*SINPH+YY*COSPH-cmpd*SINPH                                    18930
      ZSKY=-XX*CICP+yy*CISP+ZZ*SINI+CMPD*CICP                              18940
      RRHO=SQRT(YSKY**2+ZSKY**2)                                           18950
      IF(RRHO.GT.RMAX)GOTO 63                                              18960
      IF(RRHO.LT.RMIN)GOTO 126                                             18970
      THET=ARSIN(ZSKY/RRHO)                                                18980
      IF(YSKY.LT.0.) GOTO 192                                              18990
      THET=6.283185+THET                                                   19000
      GOTO 193                                                             19010
  192 THET=3.141593-THET                                                   19020
  193 IF(THET.GE.6.283185) THET=THET-6.283185                              19030
      RHHO=0.                                                              19040
      DO 52 N=1,KK                                                         19050
      ENNN=N-1                                                             19060
      ENTHET=ENNN*THET                                                     19070
   52 RHHO=RHHO+AA(N)*COS(ENTHET)+BB(N)*SIN(ENTHET)                        19080
      IF(RRHO.LE.RHHO) GOTO 80                                             19090
      SIGN=1.                                                              19100
      GOTO 81                                                              19110
   80 SIGN=-1.                                                             19120
   81 DRHO=ABS(RRHO-RHHO)                                                  19130
      IF((SIGN*OLSIGN).GE.0.) GOTO 60                                      19140
      SUMDR=DRHO+OLDRHO                                                    19150
      FACT=-(.5-DRHO/SUMDR)                                                19160
      IF(FACT.LT.0.) GOTO 198                                              19170
      RDIF=OLDIF                                                           19180
      GOTO 199                                                             19190
  198 RDIF=DIF                                                             19200
  199 CORR=FACT*RDIF*SIGN                                                  19210
      CORRR=CORR*V                                                         19220
      SUMJ=SUMJ+CORR                                                       19230
      SOMJ=SOMJ+CORRR                                                      19240
   60 IF(SIGN.LT.0.) GOTO 126                                              19250
   63 SUMJ=SUMJ+DIF                                                        19260
      SOMJ=SOMJ+DIFF                                                       19270
      xcm=xx-rmf
      dift=dif*delfi*delth
      ys=xcm*sinph+yy*cosph
      zs=-xcm*cicp+yy*cisp+zz*sini

ch    write(6,888) komp,phs,ys,zs,dift
      ilc=ilc+1         ! added by C. Hummel
      xlc(ilc)=dble(ys)
      ylc(ilc)=dble(zs)
      flc(ilc)=dble(dift)
      clc(ilc)=komp

  126 CONTINUE                                                             19280
      SOMJ=SOMJ*DELFI                                                      19290
      SUMJ=SUMJ*DELFI                                                      19300
      SOM=SOM+SOMJ                                                         19310
  136 SUM=SUM+SUMJ                                                         19320
      IF(TEST.LE..0625)GOTO 120                                            19330
      SOMHOT=SOM*DELTH                                                     19340
      SUMHOT=SUM*DELTH                                                     19350
      GOTO 121                                                             19360
  120 SUMKUL=SUM*DELTH                                                     19370
      SOMKUL=SOM*DELTH                                                     19380
  121 RETURN                                                               19390
      END                                                                  19400

