      subroutine sfcalc (ucell,istor,incel,iatcnt,bf,foc,
     & elemt,iat,rofats,nel,xdens,f000,lambda,iano,isymce,jobnam,icomp)
                                                                        STRU002
c         call sfcalc (ucell,istor,incel,iatcnt,bf,foc,
c     &   elem,iat,rofats,nel,xdens,f000,lambda,iano,isymce,jobnam)

c --- need ifried to avoid rubbish output

      parameter (maxref=500000, mst=20)                                     STRU003
                                                                        STRU004
      character*2 elemt(iat)
      character icomp*4
      character*60 phsfil,hklfil,sctfil,jobnam
      integer h,k,l                                                     STRU005            
      real cmfac(105,13),sfac(105),fp(105),fpp(105),rho(105),lambda           
      real best(mst,5),big(mst,5),ucell(istor,incel),foc(iat),bf(iat)
      real rofats(nel)                                                  STRU007
                                                                        STRU008
      data xkev /12.397639/                                             FPRI131
      pi=3.141592654                                                    STRU009
      twopi=6.2831853                                                   STRU010
      rad=180.0/pi                                                      STRU011
                                                                        STRU012
      fmax2=0.0                                                         STRU013
      prifkl=0.0
      
      hklfil = jobnam(1:index(jobnam,' ')-1)//'.hkl'                        PDB2039
      sctfil = jobnam(1:index(jobnam,' ')-1)//'.scf'                      PDB2265
      phsfil = jobnam(1:index(jobnam,' ')-1)//'.fin'

      open (11,file=sctfil, status='old')                               PDB2266
      open (12,file=hklfil, status='old')                               STRU015
      open (13,file=phsfil, status='unknown')                           STRU018
                                                                        STRU019
c --- best holds largest anomalous differences --                       STRU020
      do 1001 i=1,mst                                                   STRU021
         do 1001 j=1,5                                                  STRU022
            best(i,j)=0.0                                               STRU023
            big(i,j)=0.0                                                STRU024
 1001 continue                                                          STRU025
                                                                        STRU026
                                                                        STRU027
      write(*,'(a)')' Starting Structure factor calculation....'        STRU028
c --- read in scatfacts                                                 STRU029
      nscf=0                                                            STRU030
      do 0001 i=1,nel                                                   STRU031
         read(11,*,end=9001) (cmfac(i,j),j=1,13)                         STRU032
         nscf=nscf+1                                                    STRU033
 0001 continue                                                          STRU034
 9001 continue                                                          STRU035
      write(*,'(a,i3)')' Sets of Cromer-Mann coefficients read : ',nscf STRU036
                                                                        STRU037
c --- read in atoms  iz,x,y,z,f,b                                       STRU038
c      read(4,*,end=9002) adum                                           STRU039
c      natom=0                                                           STRU040
c      do 0003 i=1,maxat                                                 STRU041
c         read(4,*,end=9002) idum,idum,(at(i,j),j=1,6)                   STRU042
c         natom=natom+1                                                  STRU043
c 0003 continue                                                          STRU044
c 9002 continue                                                          STRU045
c      write(*,'(a,i8)')' Atom lines read : ',natom                      STRU046
                                                                        STRU047
c --- stol2 is (sin(theta)/lambda)**2                                   STRU048
      nref=0                                                            STRU049
      do 0004 i=1,maxref                                                STRU050
         read(12,*,end=9004) h,k,l,dum,stol                              STRU051
         stol2=stol*stol                                                STRU052
c ---    all sfac(iz) for this reflection stol2                         STRU053
         do 0005 j=1,nel                                                STRU054
            sfac(j)=0.0                                                 STRU055
 0005    continue                                                       STRU056
         do 0007 j=1,nscf                                               STRU057
            iz=int(cmfac(j,1))                                          STRU058
            do 0006 ii=1,4                                              STRU059
               arg=-(cmfac(j,ii+5)*stol2)                               STRU060
               sfac(iz)=sfac(iz)+cmfac(j,ii+1)*exp(arg)                 STRU061
 0006       continue                                                    STRU062
            sfac(iz)=sfac(iz)+cmfac(j,10)                               STRU063
            fp(iz)=cmfac(j,11)                                          STRU064
            fpp(iz)=cmfac(j,12)                                         STRU065
c           write(*,*) 'sfac(iz)=',iz,sfac(iz)                          STRU066
 0007    continue                                                       STRU067

c ---    for all atoms now ---                                          STRU069
         ahkl=0                                                         STRU070
         bhkl=0                                                         STRU071
c         write(*,*)h,k,l                                               STRU073
         do 0008 j=1,iatcnt                                              STRU074
            id=int(ucell(j,5))
            iz=ielect(elemt(id))                                             STRU075
c            print*, iz,id,iatcnt
c            write(*,'(a,3f10.5))') ' sfac=',sfac(iz),fp(iz),fpp(iz)     STRU076
            ph=twopi*(h*ucell(j,1)+k*ucell(j,2)+l*ucell(j,3))                    STRU077
c            write(*,'(a,f10.5)') ' ph  =',ph                            STRU078
            if (iano.eq.1) then                                         STRU079
               sfacr=sfac(iz)+fp(iz)                                    STRU080
               bxn=foc(iz)*exp(-bf(iz)*stol2)                           STRU081
               cosph=cos(ph)                                            STRU082
               sinph=sin(ph)                                            STRU083
               ahkl=ahkl+bxn*(sfacr*cosph-fpp(iz)*sinph)                STRU084
               bhkl=bhkl+bxn*(sfacr*sinph+fpp(iz)*cosph)                STRU085
            else                                                        STRU086
               fxbxn=sfac(iz)*foc(iz)*exp(-bf(iz)*stol2)                STRU087
               ahkl=ahkl+fxbxn*cos(ph)                                  STRU088
               bhkl=bhkl+fxbxn*sin(ph)                                  STRU089
            end if                                                      STRU090
 0008    continue                                                       STRU091
         fhkl=sqrt(ahkl*ahkl + bhkl*bhkl)                               STRU092
         fhkl2=fhkl*fhkl                                                STRU093
         arg=bhkl/ahkl                                                  STRU095
         phi=rad*atan2(bhkl,ahkl)                                       STRU096
         if(phi.lt.0) phi=phi+360.0
         sigf=0.5*sqrt(fhkl)                                            STRU097
         if(fhkl2.gt.fmax2) fmax2=fhkl2                                 STRU094
c ---    best anomalous signal                                          STRU098
         if (iano.eq.0) then                                            STRU099
            write(6,'(3I4,3f8.2)')h,k,l,fhkl,sigf,phi                   STRU100
         else if (mod(i,2).eq.0) then                                   STRU101
            andif=fhkl2-prhkl2                                          STRU102
            bndif=abs(andif)                                            STRU103
            if (andif.gt.0) then                                        STRU104
               andif=100.0*andif/fhkl2                                  STRU105
            else                                                        STRU106
               andif=abs(100.0*andif/prhkl2)                            STRU107
            end if                                                      STRU108
            do 2021 ii=1,mst                                            STRU109
               if (andif.gt.best(ii,5)) then                            STRU110
                  do 2001 iv=mst-1,ii,-1                                STRU111
                     do 2001 iii=1,5                                    STRU112
                        best(iv+1,iii)=best(iv,iii)                     STRU113
 2001             continue                                              STRU114
                  best(ii,5)=andif                                      STRU115
                  best(ii,1)=h                                          STRU116
                  best(ii,2)=k                                          STRU117
                  best(ii,3)=l                                          STRU118
                  best(ii,4)=fhkl2                                      STRU119
                  goto 2121                                             STRU120
               end if                                                   STRU121
 2021       continue                                                    STRU122
 2121       continue                                                    STRU123
            do 2022 ii=1,mst                                            STRU124
               if (bndif.gt.big(ii,5)) then                             STRU125
                  do 2011 iv=mst-1,ii,-1                                STRU126
                     do 2011 iii=1,5                                    STRU127
                        big(iv+1,iii)=big(iv,iii)                       STRU128
 2011             continue                                              STRU129
                  big(ii,5)=bndif                                       STRU130
                  big(ii,1)=h                                           STRU131
                  big(ii,2)=k                                           STRU132
                  big(ii,3)=l                                           STRU133
                  big(ii,4)=fhkl2                                       STRU134
                  goto 2111                                             STRU135
               end if                                                   STRU136
 2022       continue                                                    STRU137
 2111       continue                                                    STRU138
            write(13,'(3I4,3f8.2)')h,k,l,fhkl,sigf,phi           
         else                                                           STRU140
            write(13,'(3I4,3f8.2)')h,k,l,fhkl,sigf,phi                        STRU141
         end if                                                         STRU142
         prhkl2=fhkl2                                                   STRU143
         nref=nref+1                                                    STRU144
         if(mod(nref,100).eq.0) then                                    STRU145
         write(*,'(1x,a,i8,a,3I4,2f8.2)')'Refl',nref,' :',h,k,l,fhkl,phiSTRU146
c        write(*,'(1x,a,f10.0,a,f10.0)')                                STRU147
c     &   'Best rel. anomalous difference so far :',                    STRU148
c     &   best(1,5),' % at I = ',best(1,4)                              STRU149
c         write(*,'(1x,a,f10.0,a,f10.0)')                               STRU150
c     &   'Biggest anomalous difference so far   :',                    STRU151
c     &   big(1,5),'   at I = ',big(1,4)                                STRU152
         end if                                                         STRU153
 0004 continue                                                          STRU154
 9004 continue                                                          STRU155
                                                                        STRU156
      scale=1000.0                                                      STRU157
c --- print table of anomalous signals
      if (iano.eq.1) then
      write(*,*) ' Intensities scaled to ',scale                        STRU158
      write(*,*)'      Best relative differences     Biggest difference'STRU159
      write(*,*)   '  h    k    l       I      %I     h    k    l      ISTRU160
     &      dI      %I '                                                STRU161
      do 2223 i=1,mst                                                   STRU162
          best(i,4)=scale*best(i,4)/fmax2                               STRU163
c         best(i,5)=scale*best(i,5)/fmax2                               STRU164
         big(i,4)=scale*big(i,4)/fmax2                                  STRU165
         big(i,5)=scale*big(i,5)/fmax2                                  STRU166
         write(*,'(2(3f5.0,2f8.2,1x),f6.2)') (best(i,j),j=1,5),         STRU167
     &    (big(i,ii),ii=1,5),100*big(i,5)/big(i,4)                      STRU168
 2223 continue                                                          STRU169
      end if
c --- iano=1                                                                        STRU170
      close(11)                                                          STRU171
      close(12)                                                          STRU172
      close(13,status='keep')                                            STRU174
      end                                                               STRU175
                                                                        STRU176
