c	
C	LW7 model #342
C
	COMMON /FUN/ pr(768),ob(768)
	real d(768,2)
	external funk,pp
	call sas('lw7.dat','lw7342.lst','lw7342.fit')
	call assign(1,'lw7341.par')
	write(7,"(' *** LW7 model 342: mmc mds 3D ***')")
	read(9,"(16f7.4)")d
	call runstp(funk,192,768,-2,d,ob,pr,pp)
	stop
	end

	subroutine funk
	COMMON /STEP/ NV,NTRACE,MATRIX,CHISQ,MASK(1000),X(1000)
	COMMON /FUN/ pr(768),ob(768)
	real uar(3,16),uvr(3,16),xar(3,16),xvr(3,16),pa(16),pv(16)
	real uas(3,16),uvs(3,16),xas(3,16),xvs(3,16),pb(16)
	equivalence (x(1),xar),(x(49),xvr),(x(97),xas),(x(145),xvs)
	real vi(16,16),au(16,16),bi(16,16)
	equivalence (pr(1),au),(pr(257),vi),(pr(513),bi)
	integer r,s,d
	static uar,uvr,uas,uvs,rbi
	data ii/0/
c	ntrace=1
	if(ii.eq.1)go to 99
	read(1,"(3f7.4)")uar,uvr,uas,uvs
	do 5 s=1,16
	do 5 d=1,3
	xar(d,s)=uar(d,s)/8.+.5
	xas(d,s)=uas(d,s)/8.+.5
	xvr(d,s)=uvr(d,s)/8.+.5
5	xvs(d,s)=uvs(d,s)/8.+.5
	ii=1
99	call rand(12345)
	nc=1000
	dp=1./nc

	do 10 s=1,16
	do 10 d=1,3
	uar(d,s)=(xar(d,s)-.5)*8.
10	uas(d,s)=(xas(d,s)-.5)*8.
	do 20 s=1,16
	   do 50 r=1,16
50	    au(r,s)=0.
          do 20 n=1,nc
          uas1=uas(1,s)+rnorm(kk)
          uas2=uas(2,s)+rnorm(kk)
          uas3=uas(3,s)+rnorm(kk)
	   do 30 r=1,16
30	      pa(r)=tvn(uas1,uas2,uas3,uar(1,r),uar(2,r),uar(3,r))
	   pmn=10000.
           nrm=1
	   do 40 r=1,16
	      if(pa(r).gt.pmn)go to 40
	      pmn=pa(r)
	      nrm=r
40	      continue
20	    au(nrm,s)=au(nrm,s)+dp

	do 12 s=1,16
	do 12 d=1,3
	uvr(d,s)=(xvr(d,s)-.5)*8.
12	uvs(d,s)=(xvs(d,s)-.5)*8.
	do 22 s=1,16
	   do 52 r=1,16
52	    vi(r,s)=0.
          do 22 n=1,nc
          uvs1=uvs(1,s)+rnorm(kk)
          uvs2=uvs(2,s)+rnorm(kk)
          uvs3=uvs(3,s)+rnorm(kk)
	  do 32 r=1,16
32	    pv(r)=tvn(uvs1,uvs2,uvs3,uvr(1,r),uvr(2,r),uvr(3,r))
	   pmn=10000.
           nrm=1
	   do 42 r=1,16
	      if(pv(r).gt.pmn)go to 42
	      pmn=pv(r)
	      nrm=r
42	      continue
22	    vi(nrm,s)=vi(nrm,s)+dp

	do 23 s=1,16
	   do 53 r=1,16
53	    bi(r,s)=0.
          do 23 n=1,nc
          uas1=uas(1,s)+rnorm(kk)
          uas2=uas(2,s)+rnorm(kk)
          uas3=uas(3,s)+rnorm(kk)
          uvs1=uvs(1,s)+rnorm(kk)
          uvs2=uvs(2,s)+rnorm(kk)
          uvs3=uvs(3,s)+rnorm(kk)
	  do 33 r=1,16
33	    pb(r)=svn(uas1,uas2,uas3,uvs1,uvs2,uvs3,
&               uar(1,r),uar(2,r),uar(3,r),uvr(1,r),uvr(2,r),uvr(3,r))
	   pmn=10000.
           nrm=1
	   do 43 r=1,16
	      if(pb(r).gt.pmn)go to 43
	      pmn=pb(r)
	      nrm=r
43	      continue
23	    bi(nrm,s)=bi(nrm,s)+dp

	call cchi(ob(513),pr(513),256)
	rbi=sqrt(chisq/256.)
	call cchi(ob,pr,768)

	RETURN
	entry pp(is,rmsd)
	write(7,117)is,rmsd,rbi,uar,uvr,uas,uvs
117	format('Subject 'i2' rmsd,rbi= '2f7.4,16(/' uar = '3f7.4),16(/' uvr = '3f7.4)
&                           ,16(/' uas = '3f7.4),16(/' uvs = '3f7.4))
	ii=0
	return
	end


	function tvn(x,y,z,ux,uy,uz)
	fx=(x-ux)
	fy=(y-uy)
	fz=(z-uz)
	tvn=fx*fx+fy*fy+fz*fz
	return
	end

	function svn(x,y,z,w,a,b,ux,uy,uz,uw,ua,ub)
	fx=(x-ux)
	fy=(y-uy)
	fz=(z-uz)
	fw=(w-uw)
	fa=(a-ua)
	fb=(b-ub)
	svn=fx*fx+fy*fy+fz*fz+fw*fw+fa*fa+fb*fb
	return
	end

        function rnorm(n)
        common /rn/gset
        data iset/0/
        if(iset.ne.0) go to 1
2         v1=2.*rand(0)-1.
          v2=2.*rand(0)-1.
          r=v1*v1+v2*v2
          if(r.gt.1.)go to 2
        fac=sqrt(-2.*alog(r)/r)
        gset=v1*fac
        iset=1
        rnorm=v2*fac
        return
1       iset=0
        rnorm=gset
        return
        end