      subroutine mcplot
c
c     Plots MONTE CARLO results, as contained in array mca(maxitc,mxparm)
c
      include 'constnts.inc'
      include 'model.inc'
c
      real*4 x(maxitc),y(maxitc)
      character device*4,plotfile*15,title*43,parname*17,data*80
      character parname1*17,parname2*17
      character*1 digits(10)
      logical query,dogauss,foundnum,mcrun,smalle
      integer*4 nbin,parms(mxparm),indx(maxitc),mxlines
      parameter (nbin=10)
      parameter (mxlines=3000)
      parameter (mxcomb=mxparm*(mxparm-1)/2)
      integer*4 parmc(mxcomb)
      real*4 x1(nbin*2),y1(nbin*2)
      real*4 ichi(maxitc),fchi(maxitc),michi,mfchi
c
      data digits/'0','1','2','3','4','5','6','7','8','9'/
c
      mcrun=.false.
      do i=1,maxitc
       do j=1,mxparm
        if(mca(i,j).ne.0)then
         mcrun=.true.
         goto 25
        endif
       enddo
      enddo
   25 if(.not.mcrun)then
       write(outc,'(a,a,$)')' Please enter file name with MT.CARLO ',
     &     'results: '
       read(inc,30)data
       open(2,file=data,status='old',err=26)
       do i=1,maxitc            ! Clear MONTE CARLO results array
        do j=1,mxparm
         mca(i,j)=0.D0
        enddo
       enddo
       goto 27
   26  write(outc,*)'*** ERROR: File not found!'
       return
   27  do i=1,mxlines
        read(2,30,end=28)data
        if(index(data,'MONTE CARLO').ne.0)then  ! Here are MT CARLO results
         do j=1,mxlines                         ! skip header lines
          read(2,30,end=28)data
          if(index(data,'!').eq.0)goto 29
         enddo
   29    do j=1,mxlines
          if(index(data,'MONTE CARLO').ne.0)goto 291
          read(data,*,err=28,end=28)(mca(j,k),k=1,9),ichi(j),fchi(j)
          read(2,30,end=28)data
          if(index(data,'MONTE CARLO').ne.0)goto 28
          read(data,*,err=28,end=28)(mca(j,k),k=10,16)
          mca(j,13)=mca(j,13)-2440000.D0
          read(2,30,end=291)data
         enddo
        endif
       enddo
   28  write(outc,*)'*** ERROR: File corrupt!'
       close(2)
       return
  291  close(2)
       r=0.
       e=0.
       do j=1,maxitc
        if(mca(j,11).ne.0)then
         r=r+1.
         e=e+real(mca(j,11))
        endif
       enddo
       if(r.ne.0)then
        e=e/r
       else
        e=0.
       endif
       ir=r
       if(e.lt.0.01)then
        smalle=.true.
        do j=1,ir
         mca(j,11)=mca(j,11)*10.D0
        enddo
       else
        smalle=.false.
       endif
      endif
c
c     Count number of valid records
      mciter=0
      do i=1,maxitc
       irec=0
       do j=1,mxparm
        if(mca(i,j).ne.0)irec=1
       enddo
       if(irec.eq.1)then
        mciter=mciter+1
       else
        goto 24
       endif
      enddo
   24 if(mciter.eq.0)then
       write(outc,*)'*** ERROR: no valid records found!'
       return
      endif
c
      write(outc,*)'Please enter your choice: '
      write(outc,*)'-1- Histograms'
      write(outc,*)'-2- All parameter combinations'
      write(outc,*)'-3- System parameters only (diameters, colors)'
      write(outc,*)'-4- Initial and final agreement factors'
      write(outc,*)'-5- Capella special code'
      write(outc,'(a,$)')' >'
      read(inc,20)device
      if(device.eq.'2')goto 1000
      if(device.eq.'3')goto 2000
      if(device.eq.'4')goto 3000
      if(device.eq.'5')goto 4000
c-----Part 1------------------------------------------------------------
      write(outc,'(a,$)')
     &      ' Screen (s) or hardcopy (h,p,v) or exit (n) ? '
      read(inc,20)device
   20 format(a4)
      if(device.eq.'h'.or.device.eq.'H')then
       plotfile='plotmtc1.lpm'//'/lpm'
      else if(device.eq.'p'.or.device.eq.'P')then
       plotfile='plotmtc1.ps'//'/ps'
      else if(device.eq.'v'.or.device.eq.'V')then
       plotfile='plotmtc1.ps'//'/vps'
      else if(device.eq.'s'.or.device.eq.'S')then
       plotfile='/xs'
      else
       return
      endif
      if(query(' Plot Gaussians for comparison ? (y/n) '))then
       dogauss=.true.
      else
       dogauss=.false.
      endif
      nx=4
      ny=3
      innums=12
      parms(1)=1
      parms(2)=2
      do i=3,12
       parms(i)=i+4
      enddo
      do i=13,16
       parms(i)=0
      enddo
      if(query(' Plot selected parameters ? (y/n) '))then
       do i=1,mxparm
        parms(i)=0
       enddo
       innums=0
       write(outc,'(a,$)')' Please enter ids: '
       read(inc,30)data
   30  format(a80)
       foundnum=.false.
       do 40 i=1,80
       do j=1,10
        if(data(i:i).eq.digits(j))then
         foundnum=.true.
         goto 40
        endif
       enddo
       if(foundnum)then
        innums=innums+1
        foundnum=.false.
       endif
   40  continue
       read(data,*)(parms(i),i=1,innums)
       write(outc,'(a,$)')' Enter number of plots in x and y dir.: '
       read(inc,*)nx,ny
      endif
      call pgbegin(0,plotfile,nx,ny)
      call pgsch(2.0)
c
      do 100 iparm=1,mxparm
       do j=1,innums
        if(iparm.eq.parms(j))goto 101
       enddo
       goto 100
  101  continue
       xmin=real(mca(1,iparm))
       xmax=xmin
       do i=1,mciter
        x(i)=real(mca(i,iparm))
        xmin=min(xmin,x(i))
        xmax=max(xmax,x(i))
       enddo
       if(xmax.eq.xmin)then
        xmax=xmax+0.1*abs(xmax)
        xmin=xmin-0.1*abs(xmin)
        if(xmax.eq.0.and.xmin.eq.0)then
         xmax= 0.1
         xmin=-0.1
        endif
       endif
       n=mciter
       if(iparm.lt.3.or.iparm.gt.6)then  ! Forget about DR's for the moment
        call pgsci(5)
        call pghist(n,x,xmin,xmax,nbin,0)
        call pgsch(2.5)
        call pgsci(7)
        goto(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16)iparm
    1   write(parname,63)'diameter 1       '
        xpos=real(diam(1,1))
        goto 60
    2   write(parname,63)'diameter 2       '
        xpos=real(diam(2,1))
        goto 60
    3   write(parname,63)'d.ratio 12       '
        xpos=real(dratio(1,2))
        goto 60
    4   write(parname,63)'d.ratio 13       '
        xpos=real(dratio(1,3))
        goto 60
    5   write(parname,63)'d.ratio 22       '
        xpos=real(dratio(2,2))
        goto 60
    6   write(parname,63)'d.ratio 23       '
        xpos=real(dratio(2,3))
        goto 60
    7   write(parname,63)'mag.diff. (red)  '
        xpos=real(delmag(2,1))
        goto 60
    8   write(parname,63)'mag.diff. (green)'
        xpos=real(delmag(2,2))
        goto 60
    9   write(parname,63)'mag.diff. (blue) '
        xpos=real(delmag(2,3))
        goto 60
   10   write(parname,63)'maj.axis         '
        xpos=real(a_true)
        goto 60
   11   if(smalle)then
         write(parname,63)'eccentricity *10 '
         xpos=real(e_true)*10.
        else
         write(parname,63)'eccentricity     '
         xpos=real(e_true)    
        endif
        goto 60
   12   write(parname,63)'period           '
        xpos=real(period)
        goto 60
   13   write(parname,63)'epoch-2440000    '
        xpos=real(epoch)
        goto 60
   14   write(parname,63)'inclination      '
        xpos=real(inclin)
        goto 60
   15   write(parname,63)'arg.asc_node     '
        xpos=real(asc_node)
        goto 60
   16   write(parname,63)'arg.periastron   '
        xpos=real(arg_peri)
        goto 60
   60   call pgpoint(1,xpos,0.0,13)
        call pgsci(5)
        if(dogauss.and.parmerr(iparm).ne.0)then
         xinc=(xmax-xmin)/nbin
         do i=1,nbin*2
          x1(i)=xmin+(i-1)*xinc/2.
          y1(i)=1./sqrt(2.0*3.1415)
     &         *exp(-0.5*((x1(i)-xpos)/real(parmerr(iparm)))**2)
     &         *real(n)*(xinc/real(parmerr(iparm)))
         enddo
         call pgline(nbin*2,x1,y1)
        endif
        if(iparm.eq.2.and.innums.eq.0)then
         write(title,61)starno
   61    format('MT CARLO error simulation for star FK5 ',i4)
         call pgmtext('T',3.0,0.5,0.5,title)
        endif
        call pgmtext('B',2.0,0.5,0.5,parname)
c.......Now determine interval containing 2/3 of the data points
        call indexx(n,x,indx)
        dist=abs(x(indx(1))-xpos)
        imin=1
        do i=1,n
         if(abs(x(indx(i))-xpos).lt.dist)then
          imin=i
          dist=abs(x(indx(i))-xpos)
         endif
        enddo
        xlo=x(indx(max(imin-n/3,1)))
        xhi=x(indx(min(imin+n/3,n)))
        call pgsci(7)
        call pgpoint(1,xlo,0.0,28)
        call pgpoint(1,xhi,0.0,29)
        call pgsch(2.0)
        write(title,62)xlo,xhi
   62   format('2/3 int.: ',f8.3,' - ',f8.3)
        call pgsci(6)
        call pgmtext('T',1.0,0.5,0.5,title)
       endif
  100 continue

      call pgend
      return
c
                                       
   63 format(a17)                         
c-----Part 2------------------------------------------------------------
 1000 write(outc,'(a,$)')
     &      ' Screen (s) or hardcopy (h,p,v) or exit (n) ? '
      read(inc,20)device
      if(device.eq.'h'.or.device.eq.'H')then
       plotfile='plotmtc2.lpm'//'/lpm'
      else if(device.eq.'p'.or.device.eq.'P')then
       plotfile='plotmtc2.ps'//'/ps'
      else if(device.eq.'v'.or.device.eq.'V')then
       plotfile='plotmtc2.ps'//'/vps'
      else if(device.eq.'s'.or.device.eq.'S')then
       plotfile='/xs'
      else
       return
      endif
c
      nx=10
      ny=7
      i=1
      do iparm=1,mxparm
       do jparm=1,iparm-1
        parmc(i)=iparm*100+jparm
        i=i+1
       enddo
      enddo
      innums=mxcomb
      if(query(' Plot selected parameter combinations ? (y/n) '))then
       do i=1,mxcomb
        parmc(i)=0
       enddo
       innums=0
       write(outc,'(a,$)')' Please enter ids (100*id1+id2): '
       read(inc,30)data
       foundnum=.false.
       do 1040 i=1,80
       do j=1,10
        if(data(i:i).eq.digits(j))then
         foundnum=.true.
         goto 1040
        endif
       enddo
       if(foundnum)then
        innums=innums+1
        foundnum=.false.
       endif
 1040  continue
       read(data,*)(parmc(i),i=1,innums)
       do j=1,mxcomb
        jparm=mod(parmc(j),100)
        iparm=(parmc(j)-jparm)/100
        if(iparm.lt.jparm.and.parmc(j).gt.100)parmc(j)=100*jparm+iparm
       enddo
      endif
      write(outc,'(a,$)')' Enter number of plots in x and y dir.,'//
     $                    ' 0 0 for default: '
      read(inc,*)nx,ny
      if(nx.eq.0.or.ny.eq.0)then
       nx=10
       ny=7
      endif

      call pgbegin(0,plotfile,nx,ny)
c
      do 1001 iparm=1,mxparm
       do 1002 jparm=1,iparm-1
        if((iparm.le.2.or.iparm.ge.7).and.(jparm.le.2.or.jparm.ge.7))
     &                                                    then
         jcomb=100*iparm+jparm
         do j=1,innums
          if(parmc(j).eq.jcomb.or.
     &       iparm.eq.parmc(j).or.jparm.eq.parmc(j))goto 1003
         enddo 
         goto 1002
 1003    xmin=real(mca(1,iparm))
         xmax=xmin
         ymin=real(mca(1,jparm))
         ymax=ymin
         do i=1,mciter
          x(i)=real(mca(i,iparm))
          xmin=min(xmin,x(i))
          xmax=max(xmax,x(i))
          y(i)=real(mca(i,jparm))
          ymin=min(ymin,y(i))
          ymax=max(ymax,y(i))
         enddo
         if(xmax.eq.xmin)then
          xmax=xmax+0.1*abs(xmax)
          xmin=xmin-0.1*abs(xmin)
          if(xmax.eq.0.and.xmin.eq.0)then
           xmax= 0.1
           xmin=-0.1
          endif
         endif
         n=mciter
         call pgsch(2.0)
         call pgsci(5)
         call pgenv(xmin,xmax,ymin,ymax,0,0)
         goto(1101,1102,1103,1104,1105,1106,1107,1108,1109,1110,1111,
     &        1112,1113,1114,1115,1116)iparm
 1101    write(parname1,63)'diameter 1       '
         xpos=real(diam(1,1))
         goto 70
 1102    write(parname1,63)'diameter 2       '
         xpos=real(diam(2,1))
         goto 70
 1103    write(parname1,63)'d.ratio 12       '
         xpos=real(dratio(1,2))
         goto 70
 1104    write(parname1,63)'d.ratio 13       '
         xpos=real(dratio(1,3))
         goto 70
 1105    write(parname1,63)'d.ratio 22       '
         xpos=real(dratio(2,2))
         goto 70
 1106    write(parname1,63)'d.ratio 23       '
         xpos=real(dratio(2,3))
         goto 70
 1107    write(parname1,63)'mag.diff. (red)  '
         xpos=real(delmag(2,1))
         goto 70
 1108    write(parname1,63)'mag.diff. (green)'
         xpos=real(delmag(2,2))
         goto 70
 1109    write(parname1,63)'mag.diff. (blue) '
         xpos=real(delmag(2,3))
         goto 70
 1110    write(parname1,63)'maj.axis         '
         xpos=real(a_true)
         goto 70
 1111    if(smalle)then
          write(parname1,63)'eccentricity *10 '
          xpos=real(e_true)*10.
         else
          write(parname1,63)'eccentricity     '
          xpos=real(e_true)    
         endif
         goto 70
 1112    write(parname1,63)'period           '
         xpos=real(period)
         goto 70
 1113    write(parname1,63)'epoch-2440000    '
         xpos=real(epoch)
         goto 70
 1114    write(parname1,63)'inclination      '
         xpos=real(inclin)
         goto 70
 1115    write(parname1,63)'arg.asc_node     '
         xpos=real(asc_node)
         goto 70
 1116    write(parname1,63)'arg.periastron   '
         xpos=real(arg_peri)
         goto 70
   70    goto(1201,1202,1203,1204,1205,1206,1207,1208,1209,1210,1211,
     &        1212,1213,1214,1215,1216)jparm
 1201    write(parname2,63)'diameter 1       '
         ypos=real(diam(1,1))
         goto 71
 1202    write(parname2,63)'diameter 2       '
         ypos=real(diam(2,1))
         goto 71
 1203    write(parname2,63)'d.ratio 12       '
         ypos=real(dratio(1,2))
         goto 71
 1204    write(parname2,63)'d.ratio 13       '
         ypos=real(dratio(1,3))
         goto 71
 1205    write(parname2,63)'d.ratio 22       '
         ypos=real(dratio(2,2))
         goto 71
 1206    write(parname2,63)'d.ratio 23       '
         ypos=real(dratio(2,3))
         goto 71
 1207    write(parname2,63)'mag.diff. (red)  '
         ypos=real(delmag(2,1))
         goto 71
 1208    write(parname2,63)'mag.diff. (green)'
         ypos=real(delmag(2,2))             
         goto 71
 1209    write(parname2,63)'mag.diff. (blue) '
         ypos=real(delmag(2,3))
         goto 71
 1210    write(parname2,63)'maj.axis         '
         ypos=real(a_true)                 
         goto 71
 1211    if(smalle)then
          write(parname2,63)'eccentricity *10 '
          ypos=real(e_true)*10.
         else
          write(parname2,63)'eccentricity     '
          ypos=real(e_true)
         endif
         goto 71
 1212    write(parname2,63)'period           '
         ypos=real(period)
         goto 71
 1213    write(parname2,63)'epoch-2440000    '
         ypos=real(epoch)
         goto 71
 1214    write(parname2,63)'inclination      '
         ypos=real(inclin)
         goto 71
 1215    write(parname2,63)'arg.asc_node     '
         ypos=real(asc_node)
         goto 71
 1216    write(parname2,63)'arg.periastron   '
         ypos=real(arg_peri)
         goto 71
   71    call pglabel(parname1,parname2,' ')
         call pgpoint(n,x,y,-1)
         call pgsci(7)
         call pgpoint(1,xpos,ypos,17)
        endif
 1002  continue
 1001 continue
      call pgend
      return
c
c-----Part 3------------------------------------------------------------
 2000 write(outc,'(a,$)')
     &      ' Screen (s) or hardcopy (h,p,v) or exit (n) ? '
      read(inc,20)device
      if(device.eq.'h'.or.device.eq.'H')then
       plotfile='plotmtc3.lpm'//'/lpm'
      else if(device.eq.'p'.or.device.eq.'P')then
       plotfile='plotmtc3.ps'//'/ps'
      else if(device.eq.'v'.or.device.eq.'V')then
       plotfile='plotmtc3.ps'//'/vps'
      else if(device.eq.'s'.or.device.eq.'S')then
       plotfile='/xs'
      else
       return
      endif

      call pgbegin(0,plotfile,4,3)
c
      do 2001 iparm=1,9
       do 2002 jparm=1,iparm-1
        if((iparm.le.2.or.iparm.ge.7).and.(jparm.le.2.or.jparm.ge.7))
     &                                                    then
         xmin=real(mca(1,iparm))
         if(iparm.eq.8)xmin=real(mca(1,iparm)-mca(1,iparm-1))
         if(iparm.eq.9)xmin=real(mca(1,iparm)-mca(1,iparm-2))
         xmax=xmin
         ymin=real(mca(1,jparm))
         if(jparm.eq.8)ymin=real(mca(1,jparm)-mca(1,jparm-1))
         if(jparm.eq.9)ymin=real(mca(1,jparm)-mca(1,jparm-2))
         ymax=ymin
         do i=1,mciter
          x(i)=real(mca(i,iparm))
          if(iparm.eq.8)x(i)=real(mca(i,iparm)-mca(i,iparm-1))
          if(iparm.eq.9)x(i)=real(mca(i,iparm)-mca(i,iparm-2))
          xmin=min(xmin,x(i))
          xmax=max(xmax,x(i))
          y(i)=real(mca(i,jparm))
          if(jparm.eq.8)y(i)=real(mca(i,jparm)-mca(i,jparm-1))
          if(jparm.eq.9)y(i)=real(mca(i,jparm)-mca(i,jparm-2))
          ymin=min(ymin,y(i))
          ymax=max(ymax,y(i))
         enddo
         if(xmax.eq.xmin)then
          xmax=xmax+0.1*abs(xmax)
          xmin=xmin-0.1*abs(xmin)
          if(xmax.eq.0.and.xmin.eq.0)then
           xmax= 0.1
           xmin=-0.1
          endif
         endif
         n=mciter
         call pgsch(2.0)
         call pgsci(5)
         call pgenv(xmin,xmax,ymin,ymax,0,0)
         goto(2101,2102,2103,2104,2105,2106,2107,2108,2109,2110,2111,
     &        2112,2113,2114,2115,2116)iparm
 2101    write(parname1,63)'diameter 1       '
         xpos=real(diam(1,1))
         goto 80
 2102    write(parname1,63)'diameter 2       '
         xpos=real(diam(2,1))
         goto 80
 2103    write(parname1,63)'d.ratio 12       '
         xpos=real(dratio(1,2))
         goto 80
 2104    write(parname1,63)'d.ratio 13       '
         xpos=real(dratio(1,3))
         goto 80
 2105    write(parname1,63)'d.ratio 22       '
         xpos=real(dratio(2,2))
         goto 80
 2106    write(parname1,63)'d.ratio 23       '
         xpos=real(dratio(2,3))
         goto 80
 2107    write(parname1,63)'mag.diff. (red)  '
         xpos=real(delmag(2,1))
         goto 80
 2108    write(parname1,63)'color (green-red)'
         xpos=real(delmag(2,2)-delmag(2,1))
         goto 80
 2109    write(parname1,63)'color (blue-red) '
         xpos=real(delmag(2,3)-delmag(2,1))
         goto 80
 2110    write(parname1,63)'maj.axis         '
         xpos=real(a_true)
         goto 80
 2111    if(smalle)then
          write(parname1,63)'eccentricity *10 '
          xpos=real(e_true)*10.
         else
          write(parname1,63)'eccentricity     '
          xpos=real(e_true)
         endif
         goto 80
 2112    write(parname1,63)'period           '
         xpos=real(period)                 
         goto 80
 2113    write(parname1,63)'epoch-2440000    '
         xpos=real(epoch)
         goto 80
 2114    write(parname1,63)'inclination      '
         xpos=real(inclin)
         goto 80
 2115    write(parname1,63)'arg.asc_node     '
         xpos=real(asc_node)
         goto 80
 2116    write(parname1,63)'arg.periastron   '
         xpos=real(arg_peri)
         goto 80
   80    goto(2201,2202,2203,2204,2205,2206,2207,2208,2209,2210,2211,
     &        2212,2213,2214,2215,2216)jparm
 2201    write(parname2,63)'diameter 1       '
         ypos=real(diam(1,1))
         goto 81
 2202    write(parname2,63)'diameter 2       '
         ypos=real(diam(2,1))
         goto 81
 2203    write(parname2,63)'d.ratio 12       '
         ypos=real(dratio(1,2))
         goto 81
 2204    write(parname2,63)'d.ratio 13       '
         ypos=real(dratio(1,3))
         goto 81
 2205    write(parname2,63)'d.ratio 22       '
         ypos=real(dratio(2,2))
         goto 81
 2206    write(parname2,63)'d.ratio 23       '
         ypos=real(dratio(2,3))
         goto 81
 2207    write(parname2,63)'mag.diff. (red)  '
         ypos=real(delmag(2,1))
         goto 81
 2208    write(parname2,63)'color (green-red)'
         ypos=real(delmag(2,2)-delmag(2,1))
         goto 81
 2209    write(parname2,63)'color (blue-red) '
         ypos=real(delmag(2,3)-delmag(2,1))
         goto 81
 2210    write(parname2,63)'maj.axis         '
         ypos=real(a_true)
         goto 81
 2211    if(smalle)then
          write(parname2,63)'eccentricity *10 '
          ypos=real(e_true)*10.
         else
          write(parname2,63)'eccentricity     '
          ypos=real(e_true)
         endif
         goto 81
 2212    write(parname2,63)'period           '
         ypos=real(period)
         goto 81
 2213    write(parname2,63)'epoch-2440000    '
         ypos=real(epoch)
         goto 81
 2214    write(parname2,63)'inclination      '
         ypos=real(inclin)
         goto 81
 2215    write(parname2,63)'arg.asc_node     '
         ypos=real(asc_node)
         goto 81
 2216    write(parname2,63)'arg.periastron   '
         ypos=real(arg_peri)
         goto 81
   81    call pglabel(parname1,parname2,' ')
         call pgpoint(n,x,y,-1)
         call pgsci(7)
         call pgpoint(1,xpos,ypos,17)
        endif
 2002  continue
 2001 continue
      call pgend
      return
c
c-----Part 4------------------------------------------------------------
 3000 if(mcrun)then
       write(outc,*)
     & 'Sorry, the requested variables have not been initialized!'
       write(outc,*)'Please exit BINAFIT and restart!'
       return
      endif
      write(outc,'(a,$)')
     &      ' Screen (s) or hardcopy (h,p,v) or exit (n) ? '
      read(inc,20)device
      if(device.eq.'h'.or.device.eq.'H')then
       plotfile='plotmtc4.lpm'//'/lpm'
      else if(device.eq.'p'.or.device.eq.'P')then
       plotfile='plotmtc4.ps'//'/ps'
      else if(device.eq.'v'.or.device.eq.'V')then
       plotfile='plotmtc4.ps'//'/vps'
      else if(device.eq.'s'.or.device.eq.'S')then
       plotfile='/xs'
      else
       return
      endif
      call pgbegin(0,plotfile,1,1)
c
      michi=0.
      mfchi=0.
      j=1
      ymax=ichi(1)
      ymin=fchi(1)
      do while(ichi(j).ne.0)
       x(j)=real(j)
       michi=michi+ichi(j)
       mfchi=mfchi+fchi(j)
       ymax=max(ymax,ichi(j))
       ymin=min(ymin,fchi(j))
       j=j+1
      enddo
c
      j=j-1
      call pgenv(0.,x(j),ymin,ymax,0,0)
      michi=michi/real(j)
      mfchi=mfchi/real(j)
      write(title,3001)michi
 3001 format('MT CARLO: mean initial |gx|u2|d = ',f6.0)
      call pgmtext('T',3.0,0.5,0.5,title)
      write(title,3002)mfchi
 3002 format('MT CARLO: mean final   |gx|u2|d = ',f6.0)
      call pgmtext('T',1.5,0.5,0.5,title)
      call pglabel('MT CARLO iteration',' ',' ')
      call pgpoint(j,x,ichi,1)
      call pgpoint(j,x,fchi,0)
      call pgend
      return
c
c-----Part 5------------------------------------------------------------
 4000 continue
      i=1
      do while(mca(i,10).ne.0) ! while axis is non-zero
       dmod=0.62
       vtot=0.08
       btot=0.88
       dmv=real(mca(i,8))
       dmb=real(mca(i,9))
       r1=real(mca(i,1))
       r2=real(mca(i,2))
       rvf21=10.**(-dmv/2.5)
       rmv1=-2.5*log10(10.**(-vtot/2.5)/(1.+rvf21))
       rmv2=-2.5*log10(10.**(-rmv1/2.5)*rvf21)
       rbf21=10.**(-dmb/2.5)
       rmb1=-2.5*log10(10.**(-btot/2.5)/(1.+rbf21))
       rmb2=-2.5*log10(10.**(-rmb1/2.5)*rbf21)
       bv1=rmb1-rmv1
       bv2=rmb2-rmv2
       bc1=0.9923-2.606*bv1+2.004*bv1**2
       bc2=0.9923-2.606*bv2+2.004*bv2**2
       rmbol1=rmv1-dmod-bc1
       rmbol2=rmv2-dmod-bc2
       teff1=10.**(3.631+0.4966*bv1-0.4712*bv1**2)
       teff2=10.**(3.631+0.4966*bv2-0.4712*bv2**2)
       rl1=10.**((4.72-rmbol1)/2.5)
       rl2=10.**((4.72-rmbol2)/2.5)
       rr1=sqrt(1.13e15*rl1/teff1**4)
       rr2=sqrt(1.13e15*rl2/teff2**4)
c
       chir1=((rr1-r1)/real(parmerr(1)))**2
       chir2=((rr2-r2)/real(parmerr(2)))**2
       chidmv=((dmv-real(delmag(2,2)))/real(parmerr(8)))**2
       chidmb=((dmb-real(delmag(2,3)))/real(parmerr(9)))**2
c
      enddo
      return

      end

