      Program Structure

C     f77 structure.f slopy.f -lpgplot -lrecipes_f -lX11

      real mu0(999),mue(999),muAv(999),n(999),Re(999),butt,Mag(999)
      real gammln,b1,b2,bacc,rtbiss,func,nnn,redz(999),Hub,delta
      real Rb,Mub,alpha,beta,gam,logn(999),wksp1(99),wksp2(99)
      real linex(2),liney(2),linea(2),lineb(2),fitx(2),fity(2)
      real Rpear,Prob,zed,D,zD,probd,rs,probrs,tmpx(99),tmpy(99)
      real line1(2),line2(2),diff1(999),diff2(999),diff3(999)
      real Mu0Obs(999),ext,Dist,BR,BK,Bmag,kmag(999),RK(999),BR(999)
      real*8 OrthA,OrthB,OrthC,OrthD
      character zzg*3
      parameter (pi=3.141592654)
      common /DOG/ nnn
      external gammln

78      FORMAT(2F10.3)
        Hub=75.0  ! XXX (70 or 75)

C -----------------------------------------------------------------------  
C -----------------------------------------------------------------------  
11    format(10x,f5.2,1x,f4.2,f5.2,1x,f5.2)
      open (10,file='lit-data/BJ-1998.dat')
      read (10,*) zzg
      read (10,*) num
      do i=1,num
         read (10,11) Mag(i),n(i),re(i),mu0(i)
         n(i)=1.0/n(i)
         logn(i)=log10( n(i) )
         redz(i)=10**( (31.2-25.0)/5 )*Hub/299793.0  ! 31.2 distance modulus
C   ---------------------------------------------------------------------
          nnn=n(i)
          butt=1.9992*nnn - 0.3271    ! approx. (exact value below)
              if (nnn.gt.0.5) then
                 b1=0.6*2.*nnn
              elseif (nnn.le.0.5) then
                 b1=1.0E-4
              endif
              b2=1.2*2.*nnn
              bacc=1.0E-4
              butt=RTBISS(func,b1,b2,bacc)    ! exact
C   ---------------------------------------------------------------------
         re(i)=10**( re(i) ) * butt**(nnn) !scale-length to half-light Rad
         re(i)=log10( (100./Hub)*29.07*Re(i)*(1.-1./sqrt(1.+redz(i)))/
     :        (1.+redz(i)) )   ! arcsec to Kpc
         mue(i)=mu0(i)+2.5*butt*log10( exp(1.) ) 
         muAv(i)=mue(i)-2.5*log10( nnn*exp(butt)* 
     :           Exp(gammln(2.*nnn))/(butt**(2.0*nnn)) )            
         mag(i)=mag(i)-31.2   ! following Stiavelli et al. (2001)
C     In B&J98's paper they use 31.5 (assuming Virgo is at 20Mpc)
         enddo 
         close (unit=10)
C -----------------------------------------------------------------------            
C -----------------------------------------------------------------------            


C -----------------------------------------------------------------------            
C -----------------------------------------------------------------------            
21    format(23x,f6.2,2x,f5.2,3x,f5.2,9x,f5.2,29x,f5.2)
      open (20,file='lit-data/CCD-Vir-For.dat')
      read (20,*) zzg
      read (20,*) zzg
      read (20,*) idum,num2
      do i=num+1,(num+num2)
         read (20,21) Mag(i),re(i),mue(i),n(i),mu0(i)
         logn(i)=log10( n(i) )
C   ---------------------------------------------------------------------
          nnn=n(i)
          butt=1.9992*nnn - 0.3271    ! approx. (exact value below)
           if (nnn.gt.0.5) then
              b1=0.6*2.*nnn
           elseif (nnn.le.0.5) then
              b1=1.0E-4
           endif
           b2=1.2*2.*nnn
           bacc=1.0E-4
           butt=RTBISS(func,b1,b2,bacc)    ! exact
C   ---------------------------------------------------------------------  
           if (n(i).lt.14.0) then     ! XXX
              muAv(i)=mue(i)-2.5*log10( nnn*exp(butt)*
     :        Exp(gammln(2.*nnn))/(butt**(2.0*nnn)) )
           else
              muAv(i)=0.0
           endif
           mu0Obs(i)=mu0(i)
       mu0(i)=0.-2.5*log10( 10**(0.-mue(i)/2.5)*Exp(butt) )     !XXX
c       mu0Obs(i)=mu0(i) !+2.5*butt/log(10.)*(0.1/(10**(re(i))))**(1./n(i))   !XXX = mu(r=0.1)
           Mag(i)=Mag(i)+0.1  ! convert from 31.3 to 31.2
          enddo
          close (unit=20)
C -----------------------------------------------------------------------
C -----------------------------------------------------------------------


C -----------------------------------------------------------------------
C -----------------------------------------------------------------------
31      format(7x,f6.2,8x,f5.2,3x,f5.2,7x,f4.2,2x,f6.2,8x,f4.2,
     :         2x,f4.2,1x,f5.2,f5.2,1x,f4.2,1x,f4.2)
        open(30,file='lit-data/nuker-Faber.dat')
        read(30,*) idum,num3
        num4=idum-num3
        read (30,*) zzg
        read (30,*) zzg
        do i=(num+num2+1),(num+num2+idum)
           read (30,31) Dist,mu0(i),muAv(i),Re(i),Mag(i),BV,
     :                  Rb,Mub,alpha,beta,gam
           r50pc=0.1*(100./Dist)  ! ~ the rad" corresponding to 50 pc
           muav(i)=muAv(i)+BV ! B-band
           mu0Obs(i)=mu0(i)   ! V-band
           mu0(i)=mu0(i)+BV   ! B-band
c           mu0Obs(i)=0.0-2.5*log10(
c     :    10**(0.-Mub/2.5) * 2.0**((beta-gam)/alpha) *     ! = Mu at 0.1"
c     :    (0.1/Rb)**(0.-gam) * (1.+(0.1/Rb)**alpha)**    
c     :    ((gam-beta)/alpha) )  ! V-band
           mu0Obs(i)=0.0-2.5*log10(
     :    10**(0.-Mub/2.5) * 2.0**((beta-gam)/alpha) *     ! = Mu at 50 pc
     :    (r50pc/Rb)**(0.-gam) * (1.+(r50pc/Rb)**alpha)**    
     :    ((gam-beta)/alpha) )  ! V-band
           Mag(i)=Mag(i)+BV - 5.*log10(80.0/Hub) 
           Re(i)=Re(i)-3.0 + log10(80.0/Hub)  ! Faber used H_0=80 and pc
        enddo
        close (unit=30)
C -----------------------------------------------------------------------  
C -----------------------------------------------------------------------  


C -----------------------------------------------------------------------
C -----------------------------------------------------------------------
41      format(9x,f6.2,16x,f4.2,24x,f5.2)
        open(40,file='lit-data/Stiavelli.dat')
        read(40,*) num5
        read(40,*) zzg
        do i=(num+num2+idum+1),(num+num2+idum+num5)
           read (40,41) Mag(i),n(i),mu0(i)
           n(i)=1./n(i)
           logn(i)=log10( n(i) )
C   ---------------------------------------------------------------------
          nnn=n(i)
          butt=1.9992*nnn - 0.3271    ! approx. (exact value below)
           if (nnn.gt.0.5) then
              b1=0.6*2.*nnn
           elseif (nnn.le.0.5) then
              b1=1.0E-4
           endif
           b2=1.2*2.*nnn
           bacc=1.0E-4
           butt=RTBISS(func,b1,b2,bacc)    ! exact
C   ---------------------------------------------------------------------         
           mue(i)=mu0(i)+2.5*butt*log10( exp(1.) )
           muAv(i)=mue(i)-2.5*log10( nnn*exp(butt)*
     :             Exp(gammln(2.*nnn))/(butt**(2.0*nnn)) )
           mu0(i)=mu0(i)+0.9  ! (B-V) color term
           mue(i)=mue(i)+0.9  ! (B-V) color term
           muAv(i)=muAv(i)+0.9  ! (B-V) color term
           Mag(i)=Mag(i)+0.9  ! (B-V) color term
        enddo    
        close (unit=40)
C -----------------------------------------------------------------------
C -----------------------------------------------------------------------


C ------------------- Coma dE -------------------------------------------
C -----------------------------------------------------------------------
51      format(22x,f5.2,1x,f5.2,1x,f5.2,2x,f5.2,16x,f6.3)
        open(50,file='Coma-dE.dat')
        color=1.08 ! Johnson B - F666 (HST Red)
        read(50,*) num6
        do i=(num+num2+num3+num4+num5+1),
     :       (num+num2+num3+num4+num5+num6)
           read (50,51) mue(i),Re(i),n(i),Mag(i),Mu0Obs(i)
           logn(i)=log10( n(i) )           
C   ---------------------------------------------------------------------
          nnn=n(i)
          butt=1.9992*nnn - 0.3271    ! approx. (exact value below)
           if (nnn.gt.0.5) then
              b1=0.6*2.*nnn
           elseif (nnn.le.0.5) then
              b1=1.0E-4
           endif
           b2=1.2*2.*nnn
           bacc=1.0E-4
           butt=RTBISS(func,b1,b2,bacc)    ! exact
C   ---------------------------------------------------------------------
           mu0(i)=mue(i)-2.5*butt*log10( exp(1.) )
           muAv(i)=mue(i)-2.5*log10( nnn*exp(butt)*
     :           Exp(gammln(2.*nnn))/(butt**(2.0*nnn)) )
           Re(i)=(100./Hub)*29.07*Re(i)*(1.-1./sqrt(1.0231))/
     :        (1.0231)
           Re(i)=log10( Re(i) )
           Mag(i)=Mag(i)-25.-5.*log10(6925./Hub)
C  The following converts from F606R to Johnson B, and applies a redshift
C  dimming correction, and a dust-extinction correctio.   
           mu0(i)=mu0(i)+color-10*log10(1.0231)-0.024
           mu0Obs(i)=mu0Obs(i)+color-10*log10(1.0231)-0.024
           mue(i)=mue(i)+color-10*log10(1.0231)-0.024
           muAv(i)=muAv(i)+color-10*log10(1.0231)-0.024
           Mag(i)=Mag(i)+color-10*log10(1.0231)-0.024
        enddo
        close (unit=50)
C -----------------------------------------------------------------------  
C -----------------------------------------------------------------------  

C ------------------- Balcells' bulges ----------------------------------
C -----------------------------------------------------------------------
61      format(10x,f5.2,3x,f7.5,2x,f5.3,4x,f6.3,3x,f5.2,5x,f5.3,
     :         4x,f5.3,2(1x,f5.3),2x,f6.3,1x,f7.4,2x,f5.2)
        open(60,file='lit-data/Balcells.dat')
        read(60,*) num7
        read(60,*) zzg
        do i=(num+num2+num3+num4+num5+num6+1),
     :       (num+num2+num3+num4+num5+num6+num7)
           read (60,61) Dist,redz(i),ext,mue(i),mu0(i),re(i),n(i),
     :                  BR(i),RK(i),Bmag,BonD,Mu0Obs(i)
           logn(i)=log10( n(i) )           
C   ---------------------------------------------------------------------
          nnn=n(i)
          butt=1.9992*nnn - 0.3271    ! approx. (exact value below)
           if (nnn.gt.0.5) then
              b1=0.6*2.*nnn
           elseif (nnn.le.0.5) then
              b1=1.0E-4
           endif
           b2=1.2*2.*nnn
           bacc=1.0E-4
           butt=RTBISS(func,b1,b2,bacc)    ! exact
C   ---------------------------------------------------------------------
c           mu0(i)=mue(i)-2.5*butt*log10( exp(1.) )
           kMag(i)=Mue(i)-5.*log10(Re(i))-2.5*log10(
     : 2.*pi*nnn*exp(butt)*Exp(gammln(2.*nnn))/(butt**(2.*nnn)))
           kMag(i)=kMag(i)-25.-5.*log10(Dist)
c           Mag(i)=Mag(i)+BR(i)+RK(i) -0.23   ! B-band 
C           Mag(i)=Mag(i)-ext-10.0*log10(1.+redz(i))  ! already applied
           
           mu0(i)=mu0(i)+0.23
           Mu0Obs(i)=Mu0Obs(i)+0.23

           BonD=10**(BonD)  ! undo the logarithm
           BonD = 1./(1.+ 1./BonD)  ! This is now the B/T ratio
           Bmag=Bmag-2.5*log10(BonD)  ! This is now the Bulge magnitude
           Bmag=Bmag-25.-5.*log10(Dist)    ! based on the RC3 total value 
c           Bmag=Bmag-ext-10.0*log10(1.+redz(i)) ! already applied
           Mag(i)=Bmag

        enddo
        close (unit=60)
C -----------------------------------------------------------------------  
C -----------------------------------------------------------------------  


        call pgbegin(0,'?',1,1)
        call pgslw(2)                    !  Line-Width = 2 X normal
        call pgscf(2)                    !   Character Font = no. 2
        call pgsch(1.2)

        call pgvport(0.15,0.37,0.15,0.42)  
        call pgwindow(28.0,19.001,-0.6,1.999)
        call pgbox('bctns',0.0,0,'bctnsv',0.0,0)
        call pgmtext('b',2.2,0.5,0.5,'<\\gm>\\de\\u')
        call pgmtext('l',2.5,0.5,0.5,'log(R\\de\\u kpc)')
        call pgtext(27.2,1.5,'g)')
           call pgvport(0.37,0.59,0.15,0.42)
           call pgwindow(27.99,19.001,-0.6,1.999)
           call pgbox('bctns',0.0,0,'bctsv',0.0,0)
           call pgmtext('b',2.2,0.5,0.5,'\\gm\\de\\u')
           linex(1)=18.0
           linex(2)=30.0
           liney(1)=0.0-0.90
           liney(2)=liney(1)+1./3.*(linex(2)-linex(1))
           call pgline(2,linex,liney)   ! Kormendy rel.
           call pgtext(27.2,1.5,'h)')
        call pgvport(0.59,0.92,0.15,0.42)
        call pgwindow(26.5,13.0,-0.6,1.999)
        call pgbox('bctns',0.0,0,'bctsv',0.0,0)
        call pgmtext('b',2.2,0.5,0.5,'\\gm\\d0\\u')
        call pgtext(25.7,1.5,'i)')

        call pgvport(0.15,0.37,0.42,0.69)  
        call pgwindow(28.0,19.001,-0.65,1.2)
        call pgbox('bcts',0.0,0,'lbctnsv',0.0,0)
        call pgmtext('lv',2.5,0.5,0.5,'n')
        call pgtext(27.2,0.9,'d)')
           call pgvport(0.37,0.59,0.42,0.69)
           call pgwindow(27.99,19.001,-0.65,1.2)
           call pgbox('bcts',0.0,0,'lbctsv',0.0,0)
           call pgtext(27.2,0.9,'e)')
        call pgvport(0.59,0.92,0.42,0.69)
        call pgwindow(26.5,13.0,-0.65,1.2)
        call pgbox('bcts',0.0,0,'lbctsv',0.0,0)
c        line1(1)=13.0
c        line1(2)=27.0
c        line2(1)=log10(5.0)
c        line2(2)=log10(0.5)
        line1(1)=0.0
        line1(2)=30.0
        line2(1)=1.63 
        line2(2)=1.63 - line1(2)/14.
        call pgline(2,line1,line2)
        call pgtext(25.7,0.9,'f)')

        call pgvport(0.15,0.37,0.69,0.96)  
        call pgwindow(28.0,19.001,-11.0,-25.0)
        call pgbox('bcts',0.0,0,'bctnsv',0.0,0)
        call pgmtext('lv',2.5,0.5,0.5,'M\\dB\\u')
        call pgtext(27.2,-22.5,'a)')
           call pgvport(0.37,0.59,0.69,0.96)
           call pgwindow(27.99,19.001,-11.0,-25.0)
           call pgbox('bcts',0.0,0,'bctsv',0.0,0)
           call pgtext(27.2,-22.5,'b)')
        call pgvport(0.59,0.92,0.69,0.96)
        call pgwindow(26.5,13.0,-11.0,-25.0)
        call pgbox('bcts',0.0,0,'bctsv',0.0,0)
        linea(1)=0.0
        linea(2)=30.0
        lineb(1)=0.-29.5
        lineb(2)=lineb(1)+2./3.*(linea(2)-linea(1))
        call pgline(2,linea,lineb)
        call pgtext(25.7,-22.5,'c)')

        call pgsch(1.5)


C     ---- Lower Panel -----------------------
        do i=1,(num+num2+num3+num4+num5+num6)
           if (i.le.num) then 
              call pgvport(0.15,0.37,0.15,0.42)  
              call pgwindow(28.0,19.001,-0.6,1.999)
              call pgpoint(1,muav(i),Re(i),46)               
                 call pgvport(0.37,0.59,0.15,0.42)           
                 call pgwindow(27.99,19.001,-0.6,1.999)
                 call pgpoint(1,mue(i),Re(i),46)
              call pgvport(0.59,0.92,0.15,0.42)
              call pgwindow(26.5,13.0,-0.6,1.999)
              call pgpoint(1,mu0(i),Re(i),46)
           elseif (i.gt.num.AND.i.le.(num+num2)) then 
              call pgvport(0.15,0.37,0.15,0.42)  
              call pgwindow(28.0,19.001,-0.6,1.999)
              call pgpoint(1,muav(i),Re(i),228)
                 call pgvport(0.37,0.59,0.15,0.42)           
                 call pgwindow(27.99,19.001,-0.6,1.999)
                 call pgpoint(1,mue(i),Re(i),228)
             call pgvport(0.59,0.92,0.15,0.42)
             call pgwindow(26.5,13.0,-0.6,1.999)
             call pgpoint(1,mu0(i),Re(i),228)
           elseif (i.gt.(num+num2).AND.i.lt.(num+num2+num3)) then 
              call pgvport(0.15,0.37,0.15,0.42)
              call pgwindow(28.0,19.001,-0.6,1.999)
              call pgpoint(1,muav(i),Re(i),21)
                 call pgvport(0.37,0.59,0.15,0.42)           
                 call pgwindow(27.99,19.001,-0.6,1.999)
c                call pgpoint(1,mue(i),Re(i),21)
              call pgvport(0.59,0.92,0.15,0.42)
              call pgwindow(26.5,13.0,-0.6,1.999)
              call pgpoint(1,mu0(i),Re(i),21)
           elseif (i.gt.(num+num2+num3).AND.i.le.(num+num2+num3+num4))
     :   then 
              call pgvport(0.15,0.37,0.15,0.42)
              call pgwindow(28.0,19.001,-0.6,1.999)
              call pgpoint(1,muav(i),Re(i),17)
                 call pgvport(0.37,0.59,0.15,0.42)           
                 call pgwindow(27.99,19.001,-0.6,1.999)
c                 call pgpoint(1,mue(i),Re(i),17)
              call pgvport(0.59,0.92,0.15,0.42)
              call pgwindow(26.5,13.0,-0.6,1.999)
              call pgpoint(1,mu0(i),Re(i),17)
           elseif (i.gt.(num+num2+idum).AND.i.le.(num+num2+idum+num5))
     :   then
              call pgvport(0.15,0.37,0.15,0.42)
              call pgwindow(28.0,19.001,-0.6,1.999)
c              call pgpoint(1,muav(i),Re(i),13)
                 call pgvport(0.37,0.59,0.15,0.42)
                 call pgwindow(27.99,19.001,-0.6,1.999)
c                 call pgpoint(1,mue(i),Re(i),13)
             call pgvport(0.59,0.92,0.15,0.42)
             call pgwindow(26.5,13.0,-0.6,1.999)
c             call pgpoint(1,mu0(i),Re(i),13)
           elseif (i.gt.(num+num2+idum+num5).AND.i.le.
     :            (num+num2+idum+num5+num6))
     :   then
              call pgsci(2)
              call pgvport(0.15,0.37,0.15,0.42)
              call pgwindow(28.0,19.001,-0.6,1.999)
              call pgpoint(1,muav(i),Re(i),18)
                 call pgvport(0.37,0.59,0.15,0.42)
                 call pgwindow(27.99,19.001,-0.6,1.999)
                 call pgpoint(1,mue(i),Re(i),18)
             call pgvport(0.59,0.92,0.15,0.42)
             call pgwindow(26.5,13.0,-0.6,1.999)
             call pgpoint(1,mu0(i),Re(i),18)
             call pgsci(1)
           endif
        enddo


C     ---- Middle Panel -----------------------
        do i=1,(num+num2+num3+num4+num5+num6)
           if (i.le.num) then 
              call pgvport(0.15,0.37,0.42,0.69)  
              call pgwindow(28.0,19.001,-0.65,1.2)       
              call pgpoint(1,muav(i),logn(i),46)               
                 call pgvport(0.37,0.59,0.42,0.69)           
                 call pgwindow(27.99,19.001,-0.65,1.2)
                 call pgpoint(1,mue(i),logn(i),46)
              call pgvport(0.59,0.92,0.42,0.69)
              call pgwindow(26.5,13.0,-0.65,1.2)
              call pgpoint(1,mu0(i),logn(i),46)
           elseif (i.gt.num.AND.i.le.(num+num2)) then 
              call pgvport(0.15,0.37,0.42,0.69)  
              call pgwindow(28.0,19.001,-0.65,1.2)
              call pgpoint(1,muav(i),logn(i),228)
                 call pgvport(0.37,0.59,0.42,0.69)           
                 call pgwindow(27.99,19.001,-0.65,1.2)
                 call pgpoint(1,mue(i),logn(i),228)
             call pgvport(0.59,0.92,0.42,0.69)
             call pgwindow(26.5,13.0,-0.65,1.2)
             call pgpoint(1,mu0(i),logn(i),228)
           elseif (i.gt.(num+num2).AND.i.lt.(num+num2+num3)) then 
              call pgvport(0.15,0.37,0.42,0.69)
              call pgwindow(28.0,19.001,-0.65,1.2)
c              call pgpoint(1,muav(i),logn(i),21)
                 call pgvport(0.37,0.59,0.42,0.69)           
                 call pgwindow(27.99,19.001,-0.65,1.2)
c                 call pgpoint(1,mue(i),logn(i),21)
              call pgvport(0.59,0.92,0.42,0.69)
              call pgwindow(26.5,13.0,-0.65,1.2)
c              call pgpoint(1,mu0(i),logn(i),21)
           elseif (i.gt.(num+num2+num3).AND.i.le.(num+num2+num3+num4))
     :   then 
              call pgvport(0.15,0.37,0.42,0.69)
              call pgwindow(28.0,19.001,-0.65,1.2)
c              call pgpoint(1,muav(i),logn(i),17)
                 call pgvport(0.37,0.59,0.42,0.69)           
                 call pgwindow(27.99,19.001,-0.65,1.2)
c                 call pgpoint(1,mue(i),logn(i),17)
              call pgvport(0.59,0.92,0.42,0.69)
              call pgwindow(26.5,13.0,-0.65,1.2)
c              call pgpoint(1,mu0(i),logn(i),17)
           elseif (i.gt.(num+num2+idum).AND.i.le.(num+num2+idum+num5))
     :   then
              call pgvport(0.15,0.37,0.42,0.69)
              call pgwindow(28.0,19.001,-0.65,1.2)
              call pgpoint(1,muav(i),logn(i),13)
                 call pgvport(0.37,0.59,0.42,0.69)
                 call pgwindow(27.99,19.001,-0.65,1.2)
                 call pgpoint(1,mue(i),logn(i),13)
             call pgvport(0.59,0.92,0.42,0.69)
             call pgwindow(26.5,13.0,-0.65,1.2)
             call pgpoint(1,mu0(i),logn(i),13)
           elseif (i.gt.(num+num2+idum+num5).AND.i.le.
     :             (num+num2+idum+num5+num6))
     :   then
              call pgsci(2)
              call pgvport(0.15,0.37,0.42,0.69)
              call pgwindow(28.0,19.001,-0.65,1.2)
              call pgpoint(1,muav(i),logn(i),18)
                 call pgvport(0.37,0.59,0.42,0.69)
                 call pgwindow(27.99,19.001,-0.65,1.2)
                 call pgpoint(1,mue(i),logn(i),18)
             call pgvport(0.59,0.92,0.42,0.69)
             call pgwindow(26.5,13.0,-0.65,1.2)
             call pgpoint(1,mu0(i),logn(i),18)
             call pgsci(1)
           endif
        enddo


C     ---- Upper Panel -----------------------
        do i=1,(num+num2+num3+num4+num5+num6)
           if (i.le.num) then 
              call pgvport(0.15,0.37,0.69,0.96)  
              call pgwindow(28.0,19.001,-11.0,-25.0)
              call pgpoint(1,muav(i),Mag(i),46)               
                 call pgvport(0.37,0.59,0.69,0.96)           
                 call pgwindow(27.99,19.001,-11.0,-25.0)
                 call pgpoint(1,mue(i),Mag(i),46)
              call pgvport(0.59,0.92,0.69,0.96)
              call pgwindow(26.5,13.0,-11.0,-25.0)
              call pgpoint(1,mu0(i),Mag(i),46)
           elseif (i.gt.num.AND.i.le.(num+num2)) then 
              call pgvport(0.15,0.37,0.69,0.96)  
              call pgwindow(28.0,19.001,-11.0,-25.0)
              call pgpoint(1,muav(i),Mag(i),228)
                 call pgvport(0.37,0.59,0.69,0.96)           
                 call pgwindow(27.99,19.001,-11.0,-25.0)
                 call pgpoint(1,mue(i),Mag(i),228)
             call pgvport(0.59,0.92,0.69,0.96)
             call pgwindow(26.5,13.0,-11.0,-25.0)
             call pgpoint(1,mu0(i),Mag(i),228)
           elseif (i.gt.(num+num2).AND.i.lt.(num+num2+num3)) then 
              call pgvport(0.15,0.37,0.69,0.96)
              call pgwindow(28.0,19.001,-11.0,-25.0)
              call pgpoint(1,muav(i),Mag(i),21)
                 call pgvport(0.37,0.59,0.69,0.96)           
                 call pgwindow(27.99,19.001,-11.0,-25.0)
c                 call pgpoint(1,mue(i),Mag(i),21)
              call pgvport(0.59,0.92,0.69,0.96)
              call pgwindow(26.5,13.0,-11.0,-25.0)
              call pgpoint(1,mu0(i),Mag(i),21)
           elseif (i.gt.(num+num2+num3).AND.i.le.(num+num2+num3+num4))
     :   then 
              call pgvport(0.15,0.37,0.69,0.96)
              call pgwindow(28.0,19.001,-11.0,-25.0)
              call pgpoint(1,muav(i),Mag(i),17)
                 call pgvport(0.37,0.59,0.69,0.96)           
                 call pgwindow(27.99,19.001,-11.0,-25.0)
c                 call pgpoint(1,mue(i),Mag(i),17)
              call pgvport(0.59,0.92,0.69,0.96)
              call pgwindow(26.5,13.0,-11.0,-25.0)
              call pgpoint(1,mu0(i),Mag(i),17)
           elseif (i.gt.(num+num2+idum).AND.i.le.(num+num2+idum+num5))
     :   then
              call pgvport(0.15,0.37,0.69,0.96)
              call pgwindow(28.0,19.001,-11.0,-25.0)
              call pgpoint(1,muav(i),Mag(i),13)
                 call pgvport(0.37,0.59,0.69,0.96)
                 call pgwindow(27.99,19.001,-11.0,-25.0)
                 call pgpoint(1,mue(i),Mag(i),13)
             call pgvport(0.59,0.92,0.69,0.96)
             call pgwindow(26.5,13.0,-11.0,-25.0)
             call pgpoint(1,mu0(i),Mag(i),13)
           elseif (i.gt.(num+num2+idum+num5).AND.i.le.
     :             (num+num2+idum+num5+num6))
     :   then
              call pgsci(2)
              call pgvport(0.15,0.37,0.69,0.96)
              call pgwindow(28.0,19.001,-11.0,-25.0)
              call pgpoint(1,muav(i),Mag(i),18)
                 call pgvport(0.37,0.59,0.69,0.96)
                 call pgwindow(27.99,19.001,-11.0,-25.0)
                 call pgpoint(1,mue(i),Mag(i),18)
             call pgvport(0.59,0.92,0.69,0.96)
             call pgwindow(26.5,13.0,-11.0,-25.0)
             call pgpoint(1,mu0(i),Mag(i),18)
             call pgsci(1)
           endif
        enddo

        call pgend
C -------------------------------------------------------------------------



C   ==== UPPER PANEL ONLY ================================================
        call pgbegin(0,'?',1,1)
        call pgslw(2)                    !  Line-Width = 2 X normal
        call pgscf(2)                    !   Character Font = no. 2
        call pgsch(1.2)

        call pgvport(0.15,0.37,0.69,0.96)  
        call pgwindow(28.0,19.001,-11.0,-25.0)
        call pgbox('bncts',0.0,0,'bctnsv',0.0,0)
        call pgmtext('lv',2.5,0.5,0.5,'M\\dB\\u')
        call pgmtext('b',2.2,0.5,0.5,'<\\gm>\\de\\u')
        call pgtext(27.2,-22.5,'a)')
           call pgvport(0.37,0.59,0.69,0.96)
           call pgwindow(27.99,19.001,-11.0,-25.0)
           call pgbox('nbcts',0.0,0,'bctsv',0.0,0)
           call pgmtext('b',2.2,0.5,0.5,'\\gm\\de\\u')
           call pgtext(27.2,-22.5,'b)')
        call pgvport(0.59,0.92,0.69,0.96)
        call pgwindow(26.5,13.0,-11.0,-25.0)
        call pgbox('bctns',0.0,0,'bctsv',0.0,0)
        call pgmtext('b',2.2,0.5,0.5,'\\gm\\d0\\u')
        linea(1)=0.0
        linea(2)=30.0
        lineb(1)=0.-29.5
        lineb(2)=lineb(1)+2./3.*(linea(2)-linea(1))
        call pgline(2,linea,lineb)
        call pgtext(25.7,-22.5,'c)')

        call pgsch(1.5)

        do i=1,(num+num2+num3+num4+num5+num6)
           if (i.le.num) then 
              call pgvport(0.15,0.37,0.69,0.96)  
              call pgwindow(28.0,19.001,-11.0,-25.0)
              call pgpoint(1,muav(i),Mag(i),46)               
                 call pgvport(0.37,0.59,0.69,0.96)           
                 call pgwindow(27.99,19.001,-11.0,-25.0)
                 call pgpoint(1,mue(i),Mag(i),46)
              call pgvport(0.59,0.92,0.69,0.96)
              call pgwindow(26.5,13.0,-11.0,-25.0)
              call pgpoint(1,mu0(i),Mag(i),46)
           elseif (i.gt.num.AND.i.le.(num+num2)) then 
              call pgvport(0.15,0.37,0.69,0.96)  
              call pgwindow(28.0,19.001,-11.0,-25.0)
              call pgpoint(1,muav(i),Mag(i),228)
                 call pgvport(0.37,0.59,0.69,0.96)           
                 call pgwindow(27.99,19.001,-11.0,-25.0)
                 call pgpoint(1,mue(i),Mag(i),228)
             call pgvport(0.59,0.92,0.69,0.96)
             call pgwindow(26.5,13.0,-11.0,-25.0)
             call pgpoint(1,mu0(i),Mag(i),228)
           elseif (i.gt.(num+num2).AND.i.lt.(num+num2+num3)) then 
              call pgvport(0.15,0.37,0.69,0.96)
              call pgwindow(28.0,19.001,-11.0,-25.0)
              call pgpoint(1,muav(i),Mag(i),21)
                 call pgvport(0.37,0.59,0.69,0.96)           
                 call pgwindow(27.99,19.001,-11.0,-25.0)
c                 call pgpoint(1,mue(i),Mag(i),21)
              call pgvport(0.59,0.92,0.69,0.96)
              call pgwindow(26.5,13.0,-11.0,-25.0)
              call pgpoint(1,mu0(i),Mag(i),21)
           elseif (i.gt.(num+num2+num3).AND.i.le.(num+num2+num3+num4))
     :   then 
              call pgvport(0.15,0.37,0.69,0.96)
              call pgwindow(28.0,19.001,-11.0,-25.0)
              call pgpoint(1,muav(i),Mag(i),17)
                 call pgvport(0.37,0.59,0.69,0.96)           
                 call pgwindow(27.99,19.001,-11.0,-25.0)
c                 call pgpoint(1,mue(i),Mag(i),17)
              call pgvport(0.59,0.92,0.69,0.96)
              call pgwindow(26.5,13.0,-11.0,-25.0)
              call pgpoint(1,mu0(i),Mag(i),17)
           elseif (i.gt.(num+num2+idum).AND.i.le.(num+num2+idum+num5))
     :   then
              call pgvport(0.15,0.37,0.69,0.96)
              call pgwindow(28.0,19.001,-11.0,-25.0)
              call pgpoint(1,muav(i),Mag(i),13)
                 call pgvport(0.37,0.59,0.69,0.96)
                 call pgwindow(27.99,19.001,-11.0,-25.0)
                 call pgpoint(1,mue(i),Mag(i),13)
             call pgvport(0.59,0.92,0.69,0.96)
             call pgwindow(26.5,13.0,-11.0,-25.0)
             call pgpoint(1,mu0(i),Mag(i),13)
           elseif (i.gt.(num+num2+idum+num5).AND.i.le.
     :             (num+num2+idum+num5+num6))
     :   then
              call pgsci(2)
              call pgvport(0.15,0.37,0.69,0.96)
              call pgwindow(28.0,19.001,-11.0,-25.0)
              call pgpoint(1,muav(i),Mag(i),18)
                 call pgvport(0.37,0.59,0.69,0.96)
                 call pgwindow(27.99,19.001,-11.0,-25.0)
                 call pgpoint(1,mue(i),Mag(i),18)
             call pgvport(0.59,0.92,0.69,0.96)
             call pgwindow(26.5,13.0,-11.0,-25.0)
             call pgpoint(1,mu0(i),Mag(i),18)
             call pgsci(1)
           endif
        enddo

        call pgend
C -------------------------------------------------------------------------




C --------- Plot Mu_0 vs Mag  ---------------------------------------------
        call pgbegin(0,'?',1,1)
c        call pgslw(2)                    !  Line-Width = 2 X normal
        call pgscf(2)                    !   Character Font = no. 2
        call pgsch(1.5)
        call pgvport(0.15,0.7,0.15,0.9)
        call pgwindow(27.0,12.5,-11.5,-24.0)
        call pgbox('bctns',2.0,4,'bnctsv',2.0,4)
        call pgmtext('l',3.3,0.5,0.5,'M\\dB\\u [mag]')
        call pgmtext('b',2.8,0.5,0.5,
     :   '\\gm\\d0,B\\u [mag arcsec\\u-2\\d]')
c      linea and lineb are set earlier in the code.
        call pgline(2,linea,lineb)
        linea(1)=30.0
        linea(2)=0.0
        lineb(1)=-20.5
        lineb(2)=-20.5
        call pgsls(4)
        call pgslw(2)
        call pgline(2,linea,lineb)
        call pgslw(2)
        call pgsls(1)
        delta=0.0
        do i=1,(num+num2+num3+num4+num5+num6)
           if (i.le.num) then
              call pgpoint(1,mu0(i),Mag(i),46)
              diff1(i)=Mag(i)-( (2./3.)*mu0(i) - 29.5 )
              delta=delta+( diff1(i) )**2              
           elseif (i.gt.num.AND.i.le.(num+num2)) then
              call pgpoint(1,mu0(i),Mag(i),228)
              diff1(i)=Mag(i)-( (2./3.)*mu0(i) - 29.5 )
              delta=delta+( diff1(i) )**2              
           elseif (i.gt.(num+num2).AND.i.lt.(num+num2+num3)) then
              call pgpoint(1,mu0(i),Mag(i),21)
              diff1(i)=0.0
           elseif (i.gt.(num+num2+num3).AND.i.le.(num+num2+num3+num4))
     :     then
              call pgpoint(1,mu0(i),Mag(i),17)
              diff1(i)=0.0
           elseif (i.gt.(num+num2+idum).AND.i.le.(num+num2+idum+num5))
     :     then
              call pgsch(1.0)
              call pgpoint(1,mu0(i),Mag(i),13)
              call pgsch(1.5)
              diff1(i)=Mag(i)-( (2./3.)*mu0(i) - 29.5 )
              delta=delta+( diff1(i) )**2              
           elseif (i.gt.(num+num2+idum+num5)) then 
c     :           .AND.i.le.(num+num2+idum+num5+num6)) then 
              call pgsci(2)
              call pgpoint(1,mu0(i),Mag(i),18)  ! Coma
c              if (i.gt.(num+num2+num3+num4+num5+num6-3)) then 
c                 call pgpoint(1,mu0(i),Mag(i),24)
c              endif
              call pgsci(1)
              diff1(i)=Mag(i)-( (2./3.)*mu0(i) - 29.5 )
              delta=delta+( diff1(i) )**2              
           endif
        enddo
        call pgend
C ---------------------------------------------------------------------------
        delta=sqrt(delta/(num+num2+num3+num4+num5+num6))
        write(6,*) ' '
        write(6,*) 'scatter in Mag_B about the M_B-mu_0 relation'
        write(6,*) 'rms = ', delta
        write(6,*) ' '

C ------- Plot the residuals against log(n) and log(Re) -----------
        call pgbegin(0,'?',1,1)
        call pgslw(2)                    !  Line-Width = 2 X normal
        call pgscf(2)                    !   Character Font = no. 2
        call pgsch(1.5)

        call pgvport(0.15,0.48,0.15,0.6)
        call pgwindow(-0.7,1.2,-4.0,4.0)
        call pgbox('bctns',0.0,0,'bnctsv',0.0,0)
        call pgmtext('l',2.5,0.5,0.5,'\\gDM\\dB\\u')
        call pgmtext('b',2.2,0.5,0.5,'log(n)')
        linea(1)=-1.0
        linea(2)=2.0
        lineb(1)=0.0
        lineb(2)=0.0 
        call pgsls(2)
        call pgline(2,linea,lineb)
        call pgsls(1)
        do i=1,(num+num2+num3+num4+num5+num6)
           if (i.le.num) then
              call pgpoint(1,logn(i),diff1(i),46)
           elseif (i.gt.num.AND.i.le.(num+num2)) then
              call pgpoint(1,logn(i),diff1(i),228)
           elseif (i.gt.(num+num2).AND.i.lt.(num+num2+num3)) then
              continue
           elseif (i.gt.(num+num2+num3).AND.i.le.(num+num2+num3+num4))
     :     then
              continue 
           elseif (i.gt.(num+num2+idum).AND.i.le.(num+num2+idum+num5))
     :     then
              call pgsch(1.0)
              call pgpoint(1,logn(i),diff1(i),13)
              call pgsch(1.5)
           elseif (i.gt.(num+num2+idum+num5)) then 
c     :           .AND.i.le.(num+num2+idum+num5+num6)) then 
              call pgsci(2)
              call pgpoint(1,logn(i),diff1(i),18)  ! Coma
              call pgsci(1)
           endif
        enddo
c ------------------------------------------
        call pgvport(0.5,0.85,0.15,0.6)
        call pgwindow(-0.7,2.0,-4.0,4.0)
        call pgbox('bctns',0.0,0,'bmctsv',0.0,0)
        call pgmtext('b',2.2,0.5,0.5,'log(r\\de\\u)')
        call pgsls(2)
        call pgline(2,linea,lineb)
        call pgsls(1)
        do i=1,(num+num2+num3+num4+num5+num6)
           if (i.le.num) then
              call pgpoint(1,Re(i),diff1(i),46)
           elseif (i.gt.num.AND.i.le.(num+num2)) then
              call pgpoint(1,Re(i),diff1(i),228)
           elseif (i.gt.(num+num2).AND.i.lt.(num+num2+num3)) then
              continue
           elseif (i.gt.(num+num2+num3).AND.i.le.(num+num2+num3+num4))
     :     then
              continue
           elseif (i.gt.(num+num2+idum).AND.i.le.(num+num2+idum+num5))
     :     then
              call pgsch(1.0)
              continue
              call pgsch(1.5)
           elseif (i.gt.(num+num2+idum+num5)) then
c     :           .AND.i.le.(num+num2+idum+num5+num6)) then
              call pgsci(2)
              call pgpoint(1,Re(i),diff1(i),18)  ! Coma
              call pgsci(1)
           endif
        enddo
        call pgend
C -------------------------------------------------------------------------




C ---------------------log(n) vs mu, rad, mag for Coma --------------------
        call pgbegin(0,'?',1,1)
        call pgslw(2)                    !  Line-Width = 2 X normal
        call pgscf(2)                    !   Character Font = no. 2
        call pgsch(1.1)

        call pgvport(0.15,0.35,0.15,0.55)
        call pgwindow(22.5,17.001,-0.1,0.4)
        call pgbox('bctns',0.0,0,'lbnctsv',0.1,2)
        call pgmtext('lv',2.0,0.5,0.5,'n')
        call pgmtext('b',2.2,0.5,0.5,'\\gm\\d0,F606W\\u')
        call pgtext(22.0,0.34,'a)')
        ick=0
        do i=(num+num2+num3+num4+num5+1),
     :     (num+num2+num3+num4+num5+num6)
           call pgpoint(1,(mu0(i)-color),logn(i),18) 
c           if (i.ne.(num+num2+num3+num4+num5+3)) then  !XXX exclude Coma#82 = GMP2983
              ick=ick+1
              tmpx(ick)=mu0(i)-color
              tmpy(ick)=logn(i)
c           endif
        enddo
        write (6,*) ' '
        write (6,*) '   CORRELATION BETWEEN LOG(n) and Mu_0:'
        Rpear=0.
        Prob=0.
        zed=0.
        call pearsn(tmpx,tmpy,ick,Rpear,Prob,zed)
        write(6,*) 'Pearsons corr-coeff=',Rpear,'  Prob.=',Prob
        D=0.
        zd=0.
        probd=0.
        rs=0.
        probrs=0.
        call spear(tmpx,tmpy,ick,wksp1,wksp2,D,zd,probd,rs,probrs)
        write(6,*) 'Spearman rank-order corr-coeff=',
     :              rs,'  Prob.=',probrs
        write(6,*) ' No. = ',ick
c   -----------------------------
        open(45, file='TMP.DAT', status='unknown')
        write (45,*) ick
        do i=1,ick
           write(45,78) tmpy(i),tmpx(i)
        enddo
        close (unit=45)
        call slopy(OrthA,OrthB,OrthC,OrthD)
        OrthA4=1.0*OrthA    ! intercept
        OrthB4=1.0*OrthB    ! error on intercept
        OrthC4=1.0*OrthC    ! slope
        OrthD4=1.0*OrthD    ! error on slope
        write(6,*) 'REDUCED MAJ AXIS: (y & x),(Mu_0 & log(n) ):'
        write(6,*) 'Slope and error: ',OrthC4,OrthD4
        write(6,*) 'Intercept and error: ',OrthA4,OrthB4
        write(6,*) ' ** Obtained from the Jackknife sampling **'
c   -----------------------------------------------------------------

        call pgvport(0.35,0.55,0.15,0.55)
        call pgwindow(24.5,20.5,-0.1,0.4)
        call pgbox('bctns',1.0,2,'lbctsv',0.0,0)
        call pgmtext('b',2.2,0.5,0.5,'\\gm\\de,F606W\\u')
        call pgtext(24.1,0.34,'b)')
        ick=0
        do i=(num+num2+num3+num4+num5+1),
     :     (num+num2+num3+num4+num5+num6)
           call pgpoint(1,(mue(i)-color),logn(i),18)
           if (i.ne.(num+num2+num3+num4+num5+3)) then  !XXX exclude 82
              ick=ick+1
              tmpx(ick)=mue(i)
           endif
        enddo
        write (6,*) ' '
        write (6,*) '   CORRELATION BETWEEN LOG(n) and Mu_e:'
        Rpear=0.
        Prob=0.
        zed=0.
        call pearsn(tmpx,tmpy,ick,Rpear,Prob,zed)
        write(6,*) 'Pearsons corr-coeff=',Rpear,'  Prob.=',Prob
        D=0.
        zd=0.
        probd=0.
        rs=0.
        probrs=0.
        call spear(tmpx,tmpy,ick,wksp1,wksp2,D,zd,probd,rs,probrs)
        write(6,*) 'Spearman rank-order corr-coeff=',
     :              rs,'  Prob.=',probrs
        call pgsch(0.9)
        call pgtext(21.5,0.225,'3645')
        call pgsch(1.1)
        write(6,*) ' No. = ',ick


        call pgvport(0.55,0.75,0.15,0.55)
        call pgwindow(-0.35,0.299,-0.1,0.4)
        call pgbox('bctns',0.0,0,'lbctsv',0.0,0)
        call pgmtext('b',2.2,0.5,0.5,'log(R\\de\\u kpc)')
        call pgtext(-0.29,0.34,'c)')
        ick=0
        do i=(num+num2+num3+num4+num5+1),
     :     (num+num2+num3+num4+num5+num6)
           call pgpoint(1,Re(i),logn(i),18)
           if (i.ne.(num+num2+num3+num4+num5+3)) then  !XXX exclude 82
              ick=ick+1
              tmpx(ick)=Re(i)
           else
              write(6,*) ' excluding ',10**(logn(i))
           endif
        enddo
        write (6,*) ' '
        write (6,*) '   CORRELATION BETWEEN LOG(n) and log(r_e):'
        Rpear=0.
        Prob=0.
        zed=0.
        call pearsn(tmpx,tmpy,ick,Rpear,Prob,zed)
        write(6,*) 'Pearsons corr-coeff=',Rpear,'  Prob.=',Prob
        D=0.
        zd=0.
        probd=0.
        rs=0.
        probrs=0.
        call spear(tmpx,tmpy,ick,wksp1,wksp2,D,zd,probd,rs,probrs)
        write(6,*) 'Spearman rank-order corr-coeff=',
     :              rs,'  Prob.=',probrs
        write(6,*) ' No. = ',ick


        call pgvport(0.75,0.95,0.15,0.55)
        call pgwindow(-15.5,-18.99,-0.1,0.4)
        call pgbox('bctns',1.0,2,'lbctsv',0.0,0)
        call pgmtext('b',2.2,0.5,0.5,'M\\dF606W\\u')
        call pgtext(-15.8,0.34,'d)')
        ick=0
        do i=(num+num2+num3+num4+num5+1),
     :     (num+num2+num3+num4+num5+num6)
           call pgpoint(1,(Mag(i)-color),logn(i),18)
           if (i.ne.(num+num2+num3+num4+num5+3)) then  !XXX exclude 82
              ick=ick+1
              tmpx(ick)=Mag(i)-color
           endif
        enddo
        call pgsch(0.9)
        call pgtext(-15.7,0.253,'2983')
        write (6,*) ' '
        write (6,*) '   CORRELATION BETWEEN LOG(n) and Mag:'
        Rpear=0.
        Prob=0.
        zed=0.
        call pearsn(tmpx,tmpy,ick,Rpear,Prob,zed)
        write(6,*) 'Pearsons corr-coeff=',Rpear,'  Prob.=',Prob
        D=0.
        zd=0.
        probd=0.
        rs=0.
        probrs=0.
        call spear(tmpx,tmpy,ick,wksp1,wksp2,D,zd,probd,rs,probrs)
        write(6,*) 'Spearman rank-order corr-coeff=',
     :              rs,'  Prob.=',probrs
        write(6,*) ' No. = ',ick
c   -----------------------------
        open(45, file='TMP.DAT', status='unknown')
        write (45,*) ick
        do i=1,ick
           write(45,78) tmpy(i),tmpx(i)
        enddo
        close (unit=45)
        call slopy(OrthA,OrthB,OrthC,OrthD)
        OrthA4=1.0*OrthA    ! intercept
        OrthB4=1.0*OrthB    ! error on intercept
        OrthC4=1.0*OrthC    ! slope
        OrthD4=1.0*OrthD    ! error on slope
        write(6,*) 'REDUCED MAJ AXIS: (y & x),(Mag & log(n) ):'
        write(6,*) 'Slope and error: ',OrthC4,OrthD4
        write(6,*) 'Intercept and error: ',OrthA4,OrthB4
        write(6,*) ' ** Obtained from the Jackknife sampling **'
c   -----------------------------------------------------------------
        write(6,*) ' '
        call pgend 
C -------------------------------------------------------------------------




C ---------------- Mag vs log(n) ------------------------------------------
        call pgbegin(0,'?',1,1)
        call pgslw(2)                    !  Line-Width = 2 X normal
        call pgscf(2)                    !   Character Font = no. 2
        call pgsch(1.8)
        call pgvport(0.15,0.8,0.15,0.74)
        call pgwindow(-0.4,1.1,-11.5,-24.0)
        call pgbox('lbctns',0.0,0,'bnctsv',2.0,4)
        call pgmtext('vl',3.5,0.5,0.5,'M\\dB\\u [mag]')
        call pgmtext('b',2.5,0.5,0.5,'n')
        open(45, file='TMP.DAT', status='unknown')
        open(65, file='forPCA.dat')
        write (45,*) num+num2+num5+num6
        write (65,*) num+num2+num5+num6
        delta=0.0
        do i=1,(num+num2+num3+num4+num5+num6)
C          For slopy.f, we create the following file:
           if (i.le.num) then
              call pgpoint(1,logn(i),Mag(i),46)
              diff2(i)=Mag(i) - ( (0.-9.5)*logn(i) - 14.1 )
              delta=delta+( diff2(i) )**2              
              write(45,78) logn(i),Mag(i)
              write(65,*) logn(i),mu0(i),Mag(i)
           elseif (i.gt.num.AND.i.le.(num+num2)) then
              call pgpoint(1,logn(i),Mag(i),228)
              diff2(i)=Mag(i) - ( (0.-9.5)*logn(i) - 14.1 )
              delta=delta+( diff2(i) )**2              
              write(45,78) logn(i),Mag(i)
              write(65,*) logn(i),mu0(i),Mag(i)
           elseif (i.gt.(num+num2+idum).AND.i.le.(num+num2+idum+num5))
     :        then
              call pgsch(1.0)
              call pgpoint(1,logn(i),Mag(i),13)
              call pgsch(1.5)
              diff2(i)=Mag(i) - ( (0.-9.5)*logn(i) - 14.1 )
              delta=delta+( diff2(i) )**2              
              write(65,*) logn(i),mu0(i),Mag(i)
              write(45,78) logn(i),Mag(i)
           elseif (i.gt.(num+num2+idum+num5)) then 
c     :       .AND.i.le.(num+num2+idum+num5+num6)) then 
              call pgsci(2)
              call pgpoint(1,logn(i),Mag(i),18)  ! Coma
              call pgsci(1)
              diff2(i)=Mag(i) - ( (0.-9.5)*logn(i) - 14.1 )
              delta=delta+( diff2(i) )**2              
              write(45,78) logn(i),Mag(i)
              write(65,*) logn(i),mu0(i),Mag(i)
           endif
        enddo
        write (45,*) ' '
        write (45,*) 'From structure.f'
        write (65,*) ' '
        write (65,*) 'From ~graham/Coma/structure.f'
        close (unit=45)
        close (unit=65)
        
c        call slopy(OrthA,OrthB,OrthC,OrthD)
        OrthA4=1.0*OrthA    ! intercept
        OrthB4=1.0*OrthB    ! error on intercept
        OrthC4=1.0*OrthC    ! slope
        OrthD4=1.0*OrthD    ! error on slope
        write(6,*) ' '
        write(6,*) ' '
        write(6,*) 'REDUCED MAJ AXIS regression of log(n) & Abs.Mag.:'
        write(6,*) 'Slope and error: ',OrthC4,OrthD4
        write(6,*) 'Intercept and error: ',OrthA4,OrthB4
        write(6,*) ' ** Obtained from the Jackknife sampling **'

        fitx(1)=0.0-0.5
        fitx(2)=1.2
c        fity(1)=0.0 - 14.1 - 9.5*fitx(1)
c        fity(2)=0.0 - 14.1 - 9.5*fitx(2)

        fity(1)=0.0 - 14.3 - 9.4*fitx(1)
        fity(2)=0.0 - 14.3 - 9.4*fitx(2)

c        call pgsls(2)
        call pgline(2,fitx,fity)

        fity(1)=0.0 + OrthA4 + OrthC4*fitx(1)
        fity(2)=0.0 + OrthA4 + OrthC4*fitx(2)
        call pgsls(1)
c        call pgline(2,fitx,fity)
        call pgend
C ------------------------------------------------------------------------
        delta=sqrt(delta/(num+num2+num5+num6))
        write(6,*) ' '
        write(6,*) 'scatter in Mag_B about the M_B-log(n) relation'
        write(6,*) 'rms = ', delta
        write(6,*) ' '

C ------- Plot the residuals in the above diagram against Mu_0 -----------
        call pgbegin(0,'?',1,1)
        call pgslw(2)                    !  Line-Width = 2 X normal
        call pgscf(2)                    !   Character Font = no. 2
        call pgsch(1.5)
        call pgvport(0.15,0.48,0.15,0.6)
        call pgwindow(26.5,3.0,-6.0,6.0)
        call pgbox('bctns',0.0,0,'bnctsv',0.0,0)
        call pgmtext('l',2.5,0.5,0.5,'\\gDM\\dB\\u')
        call pgmtext('b',2.2,0.5,0.5,'\\gm\\d0\\u')
        linea(1)=0.0
        linea(2)=27.0
        lineb(1)=0.0
        lineb(2)=0.0 
        call pgsls(2)
        call pgline(2,linea,lineb)
        call pgsls(1)
        do i=1,(num+num2+num3+num4+num5+num6)
           if (i.le.num) then
              call pgpoint(1,mu0(i),diff2(i),46)
           elseif (i.gt.num.AND.i.le.(num+num2)) then
              call pgpoint(1,mu0(i),diff2(i),228)
           elseif (i.gt.(num+num2).AND.i.lt.(num+num2+num3)) then
              continue
           elseif (i.gt.(num+num2+num3).AND.i.le.(num+num2+num3+num4))
     :     then
              continue 
           elseif (i.gt.(num+num2+idum).AND.i.le.(num+num2+idum+num5))
     :     then
              call pgsch(1.0)
              call pgpoint(1,mu0(i),diff2(i),13)
              call pgsch(1.5)
           elseif (i.gt.(num+num2+idum+num5)) then 
c     :           .AND.i.le.(num+num2+idum+num5+num6)) then 
              call pgsci(2)
              call pgpoint(1,mu0(i),diff2(i),18)  ! Coma
              call pgsci(1)
           endif
        enddo
c ------------------------------------------
        call pgvport(0.5,0.85,0.15,0.6)
        call pgwindow(-0.7,2.0,-4.0,4.0)
        call pgbox('bctns',0.0,0,'bmctsv',0.0,0)
        call pgmtext('b',2.2,0.5,0.5,'log(r\\de\\u)')
        call pgsls(2)
        call pgline(2,linea,lineb)
        call pgsls(1)
        do i=1,(num+num2+num3+num4+num5+num6)
           if (i.le.num) then
              call pgpoint(1,Re(i),diff2(i),46)
           elseif (i.gt.num.AND.i.le.(num+num2)) then
              call pgpoint(1,Re(i),diff2(i),228)
           elseif (i.gt.(num+num2).AND.i.lt.(num+num2+num3)) then
              continue
           elseif (i.gt.(num+num2+num3).AND.i.le.(num+num2+num3+num4))
     :     then
              continue
           elseif (i.gt.(num+num2+idum).AND.i.le.(num+num2+idum+num5))
     :     then
              call pgsch(1.0)
              continue
              call pgsch(1.5)
           elseif (i.gt.(num+num2+idum+num5)) then
c     :           .AND.i.le.(num+num2+idum+num5+num6)) then
              call pgsci(2)
              call pgpoint(1,Re(i),diff2(i),18)  ! Coma
              call pgsci(1)
           endif
        enddo
        call pgend
C -------------------------------------------------------------------------




C ---------------- Mag vs log(Re) ------------------------------------------
        call pgbegin(0,'?',1,1)
        call pgslw(2)                    !  Line-Width = 2 X normal
        call pgscf(2)                    !   Character Font = no. 2
        call pgsch(1.5)
        call pgvport(0.15,0.8,0.15,0.8)
        call pgwindow(-0.7,2.0,-11.5,-24.0)
        call pgbox('lbctns',0.0,0,'bnctsv',0.0,0)
        call pgmtext('vl',2.5,0.5,0.5,'M\\dB\\u')
        call pgmtext('b',2.3,0.5,0.5,'R\\de\\u')
        open(45, file='TMP.DAT', status='unknown')
        write (45,*) num+num2+num5+num6
        delta=0.0
        do i=1,(num+num2+num3+num4+num5+num6)
C          For slopy.f, we create the following file:
           if (i.le.num) then
              call pgpoint(1,Re(i),Mag(i),46)
              diff3(i)=Mag(i) - ( (0.-3.75)*Re(i) - 16.5 )
c              delta=delta+( diff3(i) )**2              
              write(45,78) Re(i),Mag(i)
           elseif (i.gt.num.AND.i.le.(num+num2)) then
              call pgpoint(1,Re(i),Mag(i),228)
              diff3(i)=Mag(i) - ( (0.-3.75)*Re(i) - 16.5 )
c              delta=delta+( diff3(i) )**2              
              write(45,78) Re(i),Mag(i)
           elseif (i.gt.(num+num2).AND.i.lt.(num+num2+num3)) then
              call pgpoint(1,Re(i),Mag(i),21)
              diff3(i)=Mag(i) - ( (0.-3.75)*Re(i) - 16.5 )
c              delta=delta+( diff3(i) )**2
              write(45,78) Re(i),Mag(i)
           elseif (i.gt.(num+num2+num3).AND.i.le.(num+num2+num3+num4))
     :     then
              call pgpoint(1,Re(i),Mag(i),17)
              diff3(i)=Mag(i) - ( (0.-3.75)*Re(i) - 16.5 )
c              delta=delta+( diff3(i) )**2
              write(45,78) Re(i),Mag(i)
           elseif (i.gt.(num+num2+idum).AND.i.le.(num+num2+idum+num5))
     :        then
              continue 
           elseif (i.gt.(num+num2+idum+num5)) then 
c     :       .AND.i.le.(num+num2+idum+num5+num6)) then 
              call pgsci(2)
              call pgpoint(1,Re(i),Mag(i),18)  ! Coma
              call pgsci(1)
              diff3(i)=Mag(i) - ( (0.-3.75)*Re(i) - 16.5 )
c              delta=delta+( diff3(i) )**2              
              write(45,78) Re(i),Mag(i)
           endif
        enddo
        write (45,*) ' '
        write (45,*) 'From structure.f'
        close (unit=45)

        write (6,*) ' The code crashes here'
c        read(5,*) zzg
c        call slopy(OrthA,OrthB,OrthC,OrthD)
        OrthA4=1.0*OrthA    ! intercept
        OrthB4=1.0*OrthB    ! error on intercept
        OrthC4=1.0*OrthC    ! slope
        OrthD4=1.0*OrthD    ! error on slope
        write(6,*) ' '
        write(6,*) ' '
        write(6,*) 'REDUCED MAJ AXIS regression of Re(n) & Abs.Mag.:'
        write(6,*) 'Slope and error: ',OrthC4,OrthD4
        write(6,*) 'Intercept and error: ',OrthA4,OrthB4
        write(6,*) ' ** Obtained from the Jackknife sampling **'

        fitx(1)=0.0-0.7
        fitx(2)=2.5

        fity(1)=0.0 + OrthA4 + OrthC4*fitx(1)
        fity(2)=0.0 + OrthA4 + OrthC4*fitx(2)
        call pgsls(1)
c        call pgline(2,fitx,fity)  ! xxx 
        call pgend
C ------------------------------------------------------------------------
        delta=sqrt(delta/(num+num2+num5+num6))
        write(6,*) ' '
        write(6,*) 'scatter in Mag_B about the M_B-log(n) relation'
        write(6,*) 'rms = ', delta
        write(6,*) ' '

C ------- Plot the residuals in the above diagram against Mu_0 -----------
        call pgbegin(0,'?',1,1)
        call pgslw(2)                    !  Line-Width = 2 X normal
        call pgscf(2)                    !   Character Font = no. 2
        call pgsch(1.5)
        call pgvport(0.15,0.48,0.15,0.6)
        call pgwindow(26.5,3.0,-6.0,6.0)
        call pgbox('bctns',0.0,0,'bnctsv',0.0,0)
        call pgmtext('l',2.5,0.5,0.5,'\\gDM\\dB\\u')
        call pgmtext('b',2.2,0.5,0.5,'\\gm\\d0\\u')
        linea(1)=0.0
        linea(2)=27.0
        lineb(1)=0.0
        lineb(2)=0.0 
        call pgsls(2)
        call pgline(2,linea,lineb)
        call pgsls(1)
        do i=1,(num+num2+num3+num4+num5+num6)
           if (i.le.num) then
              call pgpoint(1,mu0(i),diff3(i),46)
           elseif (i.gt.num.AND.i.le.(num+num2)) then
              call pgpoint(1,mu0(i),diff3(i),228)
           elseif (i.gt.(num+num2).AND.i.lt.(num+num2+num3)) then
              call pgpoint(1,mu0(i),diff3(i),21)
           elseif (i.gt.(num+num2+num3).AND.i.le.(num+num2+num3+num4))
     :     then
c              call pgpoint(1,mu0(i),diff3(i),17)
           elseif (i.gt.(num+num2+idum).AND.i.le.(num+num2+idum+num5))
     :     then
              continue
           elseif (i.gt.(num+num2+idum+num5)) then 
c     :           .AND.i.le.(num+num2+idum+num5+num6)) then 
              call pgsci(2)
              call pgpoint(1,mu0(i),diff3(i),18)  ! Coma
              call pgsci(1)
           endif
        enddo
c ------------------------------------------
        call pgvport(0.5,0.85,0.15,0.6)
        call pgwindow(-0.7,1.6,-4.0,4.0)
        call pgbox('bctns',0.0,0,'bmctsv',0.0,0)
        call pgmtext('b',2.2,0.5,0.5,'log(n)')
        call pgsls(2)
        call pgline(2,linea,lineb)
        call pgsls(1)
        do i=1,(num+num2+num3+num4+num5+num6)
           if (i.le.num) then
              call pgpoint(1,logn(i),diff3(i),46)
           elseif (i.gt.num.AND.i.le.(num+num2)) then
              call pgpoint(1,logn(i),diff3(i),228)
           elseif (i.gt.(num+num2).AND.i.lt.(num+num2+num3)) then
              continue
           elseif (i.gt.(num+num2+num3).AND.i.le.(num+num2+num3+num4))
     :     then
              continue
           elseif (i.gt.(num+num2+idum).AND.i.le.(num+num2+idum+num5))
     :     then
              continue
           elseif (i.gt.(num+num2+idum+num5)
     :           .AND.i.le.(num+num2+idum+num5+num6)) then
              call pgsci(2)
              call pgpoint(1,logn(i),diff3(i),18)  ! Coma
              call pgsci(1)
           endif
        enddo
        call pgend
C ---------------------------------------------------------------------

C ------- Plot for Balcells paper mu_0 vs. Mag. -----------------------
        call pgbegin(0,'?',1,1)
        call pgslw(2)                    !  Line-Width = 2 X normal
        call pgscf(2)                    !   Character Font = no. 2
        call pgsch(1.5)
        call pgvport(0.15,0.5,0.55,0.95)
        call pgwindow(-12.0,-24.0,23.0,5.0)
        call pgbox('bcts',0.0,0,'bnctsv',0.0,0)
        call pgmtext('l',2.7,0.5,0.5,'\\gm\\d0,Ser\\u (H-band)')
        call pgtext(-13.0,8.0,'a)')
        do i=1,(num+num2+num3+num4+num5+num6+num7)
           if (i.le.num) then
              smut=Mu0(i)-3.5   !(B -> H, dE)
              call pgpoint(1,Mag(i),smut,46)    ! Helmut's dEs
           elseif (i.gt.num.AND.i.le.(num+num2)) then
              smut=Mu0(i)-4.0   !(B -> H, Ellip)
              call pgpoint(1,Mag(i),smut,228)  ! CCD93
           elseif (i.gt.(num+num2).AND.i.lt.(num+num2+num3)) then
c              call pgpoint(1,Mag(i),smut,21)
           elseif (i.gt.(num+num2+num3).AND.i.le.(num+num2+num3+num4))
     :     then
c              call pgpoint(1,Mag(i),smut,17)
           elseif (i.gt.(num+num2+idum).AND.i.le.(num+num2+idum+num5))
     :     then
              continue
           elseif (i.gt.(num+num2+idum+num5)
     :           .AND.i.le.(num+num2+idum+num5+num6)) then 
              smut=Mu0(i)-3.5  !(B -> H, dE)
              call pgsch(1.2)
              call pgpoint(1,Mag(i),smut,12)  ! Coma
              call pgsch(1.5)
           elseif (i.gt.(num+num2+idum+num5+num6)) then 
c     :     .AND.i.le.(num+num2+idum+num5+num6+num7)) then 
              smut=Mu0(i)
              call pgpoint(1,Mag(i),smut,17)  ! Balcells

              write (6,*) i,Mag(i),smut

           endif
        enddo
c ------------------------------------------
        call pgvport(0.15,0.5,0.15,0.55)
        call pgwindow(-12.0,-24.0,23.0,5.001)
        call pgbox('bctns',0.0,0,'bnctsv',0.0,0)
        call pgmtext('b',2.5,0.5,0.5,'M\\dB\\u [mag]')
        call pgmtext('l',2.7,0.5,0.5,'\\gm\\d0,obs\\u (H-band)')
        call pgtext(-13.0,8.0,'b)')
        do i=1,(num+num2+num3+num4+num5+num6+num7)
           if (i.le.num) then
c              call pgpoint(1,Mag(i),smut,46)   ! Helmuts dEs
           elseif (i.gt.num.AND.i.le.(num+num2)) then
              smut=Mu0Obs(i)-4.0  !(B -> H)
              call pgpoint(1,Mag(i),smut,228)   ! CCD93
           elseif (i.gt.(num+num2).AND.i.lt.(num+num2+num3)) then
              smut=Mu0Obs(i)-3.0  !(V -> H)
              call pgpoint(1,Mag(i),smut,7)    ! Nuker "power-law"
           elseif (i.gt.(num+num2+num3).AND.i.le.(num+num2+num3+num4))
     :     then
              smut=Mu0Obs(i)-3.0  !(V -> H)
              call pgpoint(1,Mag(i),smut,6)    ! Nuker "core"
           elseif (i.gt.(num+num2+idum).AND.i.le.(num+num2+idum+num5))
     :     then
              continue
           elseif (i.gt.(num+num2+idum+num5)
     :           .AND.i.le.(num+num2+idum+num5+num6)) then
              smut=Mu0Obs(i)-3.5  !(B -> H, dE)
              call pgsch(1.2)
              call pgpoint(1,Mag(i),smut,12)  ! Coma
              call pgsch(1.5)
           elseif (i.gt.(num+num2+idum+num5+num6)) then 
c     :     .AND.i.le.(num+num2+idum+num5+num6+num7)) then 
              smut=Mu0Obs(i)
              call pgpoint(1,Mag(i),smut,17)  ! Balcells
           endif
        enddo
        call pgend
C -------------------------------------------------------------------------

C ------- Test/check plot for Balcells paper: M_bh vs. K-Mag. -----------------------
        call pgbegin(0,'?',1,1)
        call pgslw(2)                    !  Line-Width = 2 X normal
        call pgscf(2)                    !   Character Font = no. 2
        call pgsch(1.5)
        call pgvport(0.15,0.5,0.55,0.95)
        call pgwindow(-17.0,-26.0,4.0,9.0)
        call pgbox('bcnts',0.0,0,'bnctsv',0.0,0)
        call pgmtext('l',2.7,0.5,0.5,'SMBH mass')
        call pgmtext('b',2.7,0.5,0.5,'Bulge K-mag')

        do i=1,(num+num2+num3+num4+num5+num6+num7)
           if (i.gt.(num+num2+idum+num5+num6)) then
              smut1=7.81 + 2.69*log10(n(i)/3.)
              call pgpoint(1,kMag(i),smut1,17)  ! Balcells
c              smut2=7.98 + 3.70*log10(n(i)/3.)
c     :        -3.10*log10(n(i)/3.)*log10(n(i)/3.) 
c              call pgpoint(1,kMag(i),smut2,27)  ! Balcells

c              smit1=10.91-0.19*(mu0(i)-0.23+RK(i)+BR(i))
c              call pgpoint(1,kMag(i),smit1,16)  ! Balcells
              smit2=8.13+0.25*(mu0(i)-0.23+RK(i)+BR(i)) -
     :0.016*(mu0(i)-0.23+RK(i)+BR(i))*(mu0(i)-0.23+RK(i)+BR(i))
              call pgpoint(1,kMag(i),smit2,6)  ! Balcells
                                          
           endif
        enddo
        call pgend

        END



C -------------------------------------------------
         real function func(b)
         real b,abe,nnn
         external gammp
         common /DOG/ nnn
            abe=2.*nnn
            func=( gammp(abe,b)-0.5 )
         return
         end


      REAL FUNCTION rtbiss(func,x1,x2,xacc)
      INTEGER JMAX
      REAL x1,x2,xacc,func,nnn
      EXTERNAL func
         common /DOG/ nnn
      PARAMETER (JMAX=100)
      INTEGER j
      REAL dx,f,fmid,xmid

c      write(6,*) 'b: ',2.*nnn,x2

      fmid=gammp(2.*nnn,x2)-0.5
      f=gammp(2.*nnn,x1)-0.5

C      write(6,*) 'c:  F-s: ',fmid,f

      if(f*fmid.ge.0.) pause 'root must be bracketed in rtbiss'
      if(f.lt.0.)then
        rtbiss=x1
        dx=x2-x1
      else
        rtbiss=x2
        dx=x1-x2
      endif
      do 11 j=1,JMAX
        dx=dx*.5
        xmid=rtbiss+dx
        fmid=gammp(2.*nnn,xmid)-0.5
        if(fmid.le.0.)rtbiss=xmid
        if(abs(dx).lt.xacc .or. fmid.eq.0.) return
 11                             continue
      pause 'too many bisections in rtbiss'
      END
