c                                  bmc.f
c                                  -----
      subroutine bmc(bmcdat,nw,win,iis,jjs,nclim,
c                ---               finaux
     1               xna,iyb,yna,bna, gna, zna,cna, 
c                        0/1  y ou b        z ou c                pour calib
c                               in            out   
     2               sundec,iswap,ipermu,icalib,fna,      
c                            si b               pour calib
     3               fl, repert)
c                    economie 

c en entree yna ou bna  (suivant iyb=0 ou 1)
c en sortie canaux de zna apres soustraction de dc et backgrounds
c        ou canaux de cna            "                   "       +calib
c------------------------------------------------------------------------

c        si iyb=1  noms de b a 38 caracteres,  b est un fits
c                                               (avant etape sequence,  "bmc")
c        si iyb=0  noms de y a 20 caracteres   b n'est pas un fits 
c                                               (avant etape flatfield, "ymz")
c
c        si icalib=1, le fichier c est calibre avec fna

      dimension xr(11,3,2),yr(11,3,2),coef(4,2),
     1          tabp(1024),tabq(1024),kz(512)
      integer*4 delb,ftime(4000),smooth,dayf,timef,win(8),winp
    
      integer*2 tab(400,400),tabxy(400,100),lec(1024),
     1          dc(1024,1024),tabx(1024,2,8),
     2          tabent(400,400),ptab(400,400),
     3          fl(400,60,8)

      character comm*80,bnom*20,bnomcourt*20,bnomlong*38,
     1 bmcdat*6,gflat*6,dark*6,data*6,end*9,fin*3
      character ftype*9,hdeb*6
      character*15 fname(4000)
      character*38 bna
      character*20 xna,yna,zna,fna,gna,cna
      dimension xna(26),yna(26),zna(26),gna(26),fna(26),
     1          bna(26),cna(26)
      Character*80 Repert

      nclim=1
c      nseqdir=0
      smooth=1    

        do i=1,512
        kz(i)=0
        enddo
c
c
      type *,bmcdat
      rewind(8)
      call SNAME(8,bmcdat,ier)
1     format(a)
2     format(1xa)
3     format(12i6)
4     format(2x12i6)
5     format(i7)
      read(8,1)comm
      write(3,2)comm

        do i=1,8 
      read(8,3)winp,im,jmp,n1,n2,maxlev,idc,iscat,smooth,lbda,delb,
     1         invern
      write(3,4)winp,im,jmp,n1,n2,maxlev,idc,iscat,smooth,lbda,delb,
     1         invern
        if(winp.eq.win(nw))goto10
        enddo

10      if(smooth.eq.0)smooth=1
      facmax=float(maxlev)/1000.
      if(maxlev.eq.0)facmax=1000.

12    format(48x,2i6)

      iux=23
      call openold20(xna(nw),sundec,23)
c          ---------                             x  en  23

      iug=24
      call openold20sf(gna(nw),sundec,24)
c          -----------                           g  en  24

      iuy=25
      iuz=26
      iub=25
      iuc=26

      if(iyb.eq.0)then
      call openold20(yna(nw),sundec,25)
c          ---------                             y  en  25 
      call opennew20(zna(nw),sundec,26)
c          ---------                             z  en  26

      else
      call openold38(bna(nw),sundec,25, repert)
c          ---------                                b en 25
      call opennew20(cna(nw),sundec,26)
c          ---------                                c en 26

         type*,' cna ouverts'
      endif

      if(icalib.eq.1.and.iyb.eq.1)then
      iuf=29
      call openold20(fna(nw),sundec,29)
c          ---------                             f  en  29
      endif
c---------------------------  lecture fichier g
      do n=1,6
      read(iug,14)
      enddo
c        read(iug,12)jcor1,jcor2
        jcor1=0
        jcor2=0
      do n=1,13
      read(iug,14)
      enddo
c
13    format(11f7.1)
14    format(1x)
      read(iug,13)xr  !!!! 2 -> iug
      read(iug,13)yr   !!! 2 -> iug

c      type*,' bmc: yr ',yr
c
c      if(n1.ne.0.and.n2.ne.0)goto16
c      n1=1
c      n2=11
c      if(xr(11,1,1).eq.0.)n2=10
c      if(xr(10,1,1).eq.0.)n2=9
c      if(xr(9,1,1).eq.0.)n2=8

16    nm=n2-n1+1 
c     --
c
      nc=(n1+n2)/2
      xlong=0.5*(xr(nc,3,1)-xr(nc,1,1)+xr(nc,3,2)-xr(nc,1,2))
      ylarg=0.5*(yr(nc,1,2)-yr(nc,1,1)+yr(nc,3,2)-yr(nc,3,1))
      yespace=0.5*(yr(nc+1,2,1)-yr(nc,2,2)+yr(nc,2,1)-yr(nc-1,2,2))
        if(im.ne.0.and.jmp.ne.0)then
        jm=jmp
        endif
        if(im.ne.0.and.jmp.eq.0)then
        jm=im*ylarg/xlong
        endif
        if(im.eq.0.and.jmp.ne.0)then
        jm=jmp
        im=jm*xlong/ylarg
        endif
        if(im.eq.0.and.jmp.eq.0)then
        piv=xlong/ylarg
         if(piv.ge.8.)then
          im=1024
          jm=im*ylarg/xlong
          endif
          if(piv.lt.8.)then
          jm=128
          im=jm*xlong/ylarg
          endif
        endif
      if(im.gt.1024)im=1024
      if(jm.gt.128)jm=128
      ism=im/smooth
      jsm=jm/smooth

      jespace=yespace*jm/ylarg   !!! coordonnees nouvelles de z ou c

18    format(' im,jm,ism,jsm,jespace  ='5i6)
      write(3,18)im,jm,ism,jsm,jespace
c-----------------------------------------
      do i=1,1024
      do j=1,1024
      dc(i,j)=0
      enddo
      enddo
c
c DARK CURRENT
      if(idc.ne.0)then
        rewind(iux)
        read(iux)kz
        type 4,(kz(i),i=1,8)
        is=kz(2)
        js=kz(3)
          do j=1,js
          read(iux)(dc(i,j),i=1,is)
            if((j-jcor1)*(j-jcor2).lt.0)then
            do i=1,is
            dc(i,j)=dc(i,j-1)
            enddo
            endif
          enddo
      endif
c
c
C-------------------------------
      if(iyb.eq.0)then
c        --------                       Lecture du fichier y 
       type *,'y : ',yna(nw)
       rewind(iuy)  
c                    header  fichier y
        read(iuy)kz

        is=kz(2)
        js=kz(3)
c       --
        idate=kz(6)!!!!!!!!!!
        itime=kz(7)!!!!!!!!!!!!
        read(yna(nw)(2:7),'(i6)')idate
        read(yna(nw)(9:14),'(i6)')itime
        type *,'idate : ',idate !!!!!!!!!!!
        type *,'itime',itime !!!!!!!!!!!!
c        iline=kz(8)
        read(yna(nw)(19:19),'(i1)')iline
        type *,' ymz: iline:',iline

       else
c      ----                                         lecture fichier c
c Fichier d'entree est un fichier B donc un fits a acces direct:
c                                                   header + tableau
        type *,'b : ',bna(nw)       
        read(bna(nw)(17:22),'(i6)')idate
        read(bna(nw)(24:29),'(i6)')itime
        read(bna(nw)(19:19),'(i1)')iline
c        type *,'idate : ',idate
c        type *,'itime :',itime
c        type *,'iline :',iline
c Lecture du nombre de blocs headers
        call comptehead(iub,inbhead)
        call readifits(iub,'NAXIS1  ',naxis1,1)
        call readifits(iub,'NAXIS2  ',naxis2,1)
        type *,' bmc: naxis1: ',naxis1
        type *,' bmc: naxis2 : ',naxis2
        type *,'inbhead : ',inbhead
        type *,'nclim : ',nclim

      naxis1=iis
      naxis2=jjs

      if(ipermu.eq.0)then
        call readtabfits(iub,inbhead,iswap,naxis1,naxis2,nclim,tabent)
c            -----------
      else
        call readtabfits(iub,inbhead,iswap,naxis2,naxis1,nclim,ptab)
c            -----------
c        ipiv=naxis1
c        naxis1=naxis2
c        naxis2=ipiv

        do i=1,naxis1
          ip=naxis1+1-i
        do j=1,naxis2
          tabent(ip,j)=ptab(j,i)
        enddo
        enddo        

      endif
               is=naxis1
               js=naxis2

      endif
c------------------------------
c Chargement de TAB(i,j)

      if(iyb.eq.0) then
c        --------
c cas de fichier sequentiel y
                                type*,'z:  is,js ',is,js
      do j=1,js
       read(iuy)(lec(i),i=1,is)
        do i=1,is
          if(idc.ne.0)then
          kpiv=lec(i)-dc(i,j)
c                    --------
          else
          kpiv=lec(i)
          endif
        tab(i,j)=kpiv
c       --------
        enddo
      enddo
      else
c     ----
c Cas de fichier fits direct
 
          write(3,*)' bmc: is,js ',is,js
       do j=1,js
       do i=1,is	
         if(idc.ne.0)then
         kpiv=tabent(i,j)-dc(i,j)
c                        --------
         else
         kpiv=tabent(i,j)
         endif
       tab(i,j)=kpiv
c      --------
       enddo

         jpiv=js/2
       if(j.eq.js/2)write(3,*)' bmc: tabent(i tous 50,js/2),dc,tab ',
     1    (tabent(i,jpiv),i=1,is,50),(dc(i,jpiv),i=1,is,50),
     2    (tab(i,jpiv),i=1,is,50)

       enddo
      endif
c                 ------------             
c      do j=1,js
c      do i=3,is-2
c        piv1=tab(i-2,j)*facmax
c        piv2=tab(i+2,j)*facmax
c        piv=tab(i,j)
c        if(piv.gt.piv1.and.piv.gt.piv2)tab(i,j)=0.5*(piv1+piv2)
c      enddo
c      enddo

c      if(jcor1*jcor2.ne.0)then
c        do j=jcor1+1,jcor2-1
c        alp=(j-jcor1)/(jcor2-jcor1)
c        bet=1.-alp
c          do i=1,is
c          tab(i,j)=alp*tab(i,jcor2)+bet*tab(i,jcor1)
c          enddo
c        enddo
c      endif
c-------------------------------------------
c                              lectures pour icalib=1
c                             avant d'entrer dans la boucle canaux

c                                   Lecture de fl
      if(icalib.eq.1.and.iyb.eq.1)then
      rewind(iuf)
      read(iuf)na3,imf,jmf,nmf  ,imm,jmm,nmm,ndf,dayf,timef,ilinef,
c                  -----------
     1          lbda,ncf,jt100,ja100,jb100,jz100,kdec,kminm
      type *,iuf
      type *,'apres read(iuf) lbda:',lbda!!!!!!!!!!!!!
      type *,'dayf,timef',dayf,timef
      type *,'jz100',jz100
      write(3,12)iuf,imf,jmf,nmf,dayf,timef,ilinef,lbda,
     1          ncf,jt100,ja100,jb100,jz100,kdec,kminm

      ltrjf=jt100
c     ncf est pris dans le header
c      ncf=(1+nmf)/2
c             average line-center = n=ncf and
c                                   j=(ja100+jz100)/100 in i=1
c                                            jz100/100     i=(1+imf)/2
c                                     (jb100+jz100)/100    i=imf
         do n=1,nmf
         do j=1,jmf
         read(iuf)(fl(i,j,n),i=1,imf)
         enddo
         enddo
      endif

c----------------------
c                                Ecriture header fichier z ou c

      kz(1)=8 !!!!!!!!!!!!!!!! 5
      kz(2)=im
      kz(3)=jm
      kz(4)=nm
      kz(5)=im  !! a supprimer
      kz(6)=nm  !! a supprimer
      kz(7)=2   !!  .....
      kz(8)=1
      kz(9)=1
      kz(10)=4
      kz(11)=idate
      kz(12)=itime
      kz(13)=iline
      kz(14)=lbda

      if(iyb.eq.1)then
      kz(10)=10
      kz(15)=ncf
      kz(16)=jt100
      kz(17)=ja100
      kz(18)=jb100
      kz(19)=jz100
      kz(20)=kdec
      kz(21)=kminm
      endif

      type *,'kz(1a20) :', (kz(i),i=1,21)
      write(iuc)kz
c           ---
c----------------------------------------------
c                                               boucle canaux
      do90 np=n1,n2
        n=np
        if(invern.ne.0)n=n2-np+n1
c                                              n initial non inverse
c                                              np final  inverse eventuellement
c Calcul de la matrice de transformation
c      write(3,97) !!!!!!!!!!!!!!!!!!!
c      type*,' bmc: im,jm,n ',im,jm,n
      call COEFF(im,jm,xr,yr,n,coef)
c      type *,' bmc: coef ',coef
c-------------------------
c                     Interpolation du fichier C      
        do jj=1,jm
          do ii=1,im
          call PIX(ii,jj,coef,i,j,di,dj)
          tabxy(ii,jj)=(tab(i,j)*(1.-di)+tab(i+1,j)*di)*(1.-dj)
     1            +(tab(i,j+1)*(1.-di)+tab(i+1,j+1)*di)*dj
          enddo
        enddo

      write(9,*),' bmc calcule z ou c: nw, np, tab(187,25,np) ',
     1                                 nw, np, tabxy(187,25)


c      ii=im/2
c      write(3,*)' bmc: n=',n,' tabxy(im/2,j tous 10)'
c      write(3,*)(tabxy(ii,jj),jj=1,jm,10)

c-------------------------
c                          Calcul des backgrounds 
c                                 (deux courbes par canal)

      call COEFF(im,jm,xr,yr,n,coef)
c          -----             -
        jesp2=jespace/2  !!     coordonnees nouvelles de z ou c

c      if(n.eq.4)then
        jmesp=jm+jespace
        ii=im/2
      do l=1,jmesp
      jj=l-jesp2
      call pix(ii,jj,coef,i,j,di,dj)
      tabp(l)=tab(i,j)
      enddo
c        write(3,*)' coupe a im/2 ',(tabp(l),l=1,jmesp)
c      endif

        jj=1-jesp2
      do jgd=1,2
          do ii=1,im
          call PIX(ii,jj,coef,i,j,di,dj)
c        if(ii.eq.im/2)type*,' n,jgd,ii,jj, i,j ',n,jgd,ii,jj,i,j
          tabp(ii)=tab(i,j)
          enddo

          call parafit(tabp,1,im,20,tabq)
            do ii=1,im
            tabx(ii,jgd,np)=tabq(ii)+0.5
c                       --
            enddo

        jj=jm+jesp2
      enddo

c-------------------
c                  lissage, soustract. fond, et ecriture dans le cas icalib=0
      if(icalib.eq.0.and.iyb.eq.0)then
c       call ssmooth(tabxy,1024,200,ism,jsm,smooth)  !!! remettre eventuellement
c           -------

        if(iscat.ne.0)then
      do i=1,im
        kpiv=0.5*(tabx(i,1,np)+tabx(i,2,np))
          do j=1,jm
          tabxy(i,j)=tabxy(i,j)-kpiv
          enddo
      enddo
        endif

c      ii=im/2
c      write(3,*)' bmc: n=',n,' tabxy(im/2,j) apres soustract. fond'
c      write(3,*)(tabxy(ii,jj),jj=1,jm)

      do j=1,jm
      write(iuc)(tabxy(i,j),i=1,im)
      enddo

      write(9,*)' bmc ecrit z: nw,np, tabxy(187,25) ',
     1                         nw,np, tabxy(187,25)

      endif
c---------------------------------------------
c                                          Calibration
      if(icalib.eq.1)then

c  soustraction des backgrounds et division par fl 
c      type *,'etape soustraction du fond' !!!!!!!!!!!!!!!!

        fscat=1.
        if(iscat.eq.0)fscat=0.
       do j=1,jm
          do i=1,im
          piv=tabxy(i,j)-0.5*(tabx(i,1,np)+tabx(i,2,np))*fscat
          pivd=fl(i,j,np) 

c       if(j.eq.jm/2)then
c      if(i.eq.10)type*,' bmc: i=10 tabxy,tabxy-tabx,fl ',
c     1                             tabxy(i,j),piv,pivd
c      if(i.eq.300)type*,' bmc: i=300 tabxy,tabxy-tabx,fl ',
c     1                               tabxy(i,j),piv,pivd
c       endif

          tabxy(i,j)=32000
            if(pivd.ne.0)then
            pivt=piv*10000./pivd+0.5
c             pivt=piv
              if(abs(pivt).le.32000.)tabxy(i,j)=pivt

            endif
          enddo
        enddo

      write(9,*)' bmc ecrit c: nw,np, fl(187,25,np),tabxy(187,25) ',
     1                         nw,np, fl(187,25,np),tabxy(187,25)

c                                                lissage et ecriture
c       call ssmooth(tabxy,1024,200,im,jm,smooth)
c           -------
      do j=1,jm
      write(iuc)(tabxy(i,j),i=1,im)
c      if(j.eq.jm/2)type *,tabxy(150,20)
      enddo

      endif

c--------------------     fin calib du canal

c             Il n'est pas necessaire de rajouter les backgrounds
c             cmd ne les utilisera pas.

c            Ecriture des backgrounds dans le cas icalib=0 
c          if(icalib.eq.0)then
c        write(iuc)(tabx(i,1,n),i=1,im)
c        write(iuc)(tabx(i,2,n),i=1,im)
c          endif

90    continue
c                                       fin boucle canaux
c-----------------------------------------
c      if(delb.ne.0)call system('rm '//bna(nw))  !! annule par precaution!!!!!
      close(unit=23)
      close(unit=24)
      close(unit=25)
      close(unit=26)

      if(icalib.eq.1.and.iyb.eq.1)close(unit=29)

      return
      end
c********************************************************************
      subroutine COEFF(im,jm,xr,yr,n,coef)
      dimension coef(4,2),xr(11,3,2),yr(11,3,2)
      xl=im
      yl=jm
c
      x1=xr(n,1,1)
      x2=xr(n,3,1)
      x3=xr(n,1,2)
      x4=xr(n,3,2)
      y1=yr(n,1,1)
      y2=yr(n,3,1)
      y3=yr(n,1,2)
      y4=yr(n,3,2)
c
      coef(1,1)=x1
      coef(2,1)=(x2-x1)/xl
      coef(3,1)=(x3-x1)/yl
      coef(4,1)=(x1+x4-x2-x3)/(xl*yl)
      coef(1,2)=y1
      coef(2,2)=(y2-y1)/xl
      coef(3,2)=(y3-y1)/yl
      coef(4,2)=(y1+y4-y2-y3)/(xl*yl)
c
      return
      end
c--------------
      subroutine PIX(ii,jj,coef,i,j,di,dj)
      dimension coef(4,2)
      xi=coef(1,1)+ii*coef(2,1)+jj*coef(3,1)+ii*jj*coef(4,1)-0.5
      yj=coef(1,2)+ii*coef(2,2)+jj*coef(3,2)+ii*jj*coef(4,2)-0.5
      i=xi
      j=yj
      di=xi-i
      dj=yj-j
      return
      end      
c----------------------------------------------------

       Subroutine readtabfits(iu,inbhead,iswap,is,js,nclim,tabent)
c                 -----------
       integer*2 tabent(400,400)
       integer*2 sort(1440),in(1),out(1)
c Lecture de tabent(is,js)
        n=inbhead
        i=0
        j=1
        k=1
100      n=n+1
        read(iu,rec=n,iostat=ios)  (sort(t),t=1,1440)
c        type *,'apres read n =',n
        if(ios.lt.0) go to 1000
        t=1
200     i=i+1
        if(i.le.is) go to 500
         if(j.lt.js) then
           j=j+1
           i=1
         else
          goto1000
c            if(k.lt.nclim) then
c               k=k+1
c               j=1
c               i=1
c            else
c               go to 1000
c            endif
          endif

500        tabent(i,j)= sort(t)
           if(iswap.eq.1)then
              in(1)=sort(t)
              call swap(in,1,out)
              tabent(i,j)=out(1)
           endif

       t=t+1

       if (t.eq.1441) go to 100
       go to 200
1000   type *,'lecture des donnees en I2 effectuee'

          return
          end
c----------------------------------                 
