c	
C	LC1 model #342
C
	COMMON /FUN/ pr(1452),ob(1452)
	real d(1452,2)
	external funk,pp
	call sas('lc1.dat','lc1342.lst','lc1342.fit')
	call assign(1,'lc1341.par')
	write(7,"(' *** LC1 model 342: mmc mds 3D ***')")
	read(9,"(8x,f7.4)")d
	call crn
	call runstp(funk,264,1452,-2,d,ob,pr,pp)
	stop
	end

	subroutine funk
	COMMON /STEP/ NV,NTRACE,MATRIX,CHISQ,MASK(1000),X(1000)
	COMMON /FUN/ pr(1452),ob(1452)
	real uar(3,22),uvr(3,22),xar(3,22),xvr(3,22),pa(22),pv(22)
	real uas(3,22),uvs(3,22),xas(3,22),xvs(3,22),pb(22)
	equivalence (x(1),xar),(x(67),xvr),(x(133),xas),(x(199),xvs)
	real vi(22,22),au(22,22),bi(22,22)
	equivalence (pr(1),au),(pr(485),vi),(pr(969),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,22
	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 rrn
	nc=1000
	dp=1./nc

	do 10 s=1,22
	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,22
	   do 50 r=1,22
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,22
30	      pa(r)=tvn(uas1,uas2,uas3,uar(1,r),uar(2,r),uar(3,r))
	   pmn=10000.
           nrm=1
	   do 40 r=1,22
	      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,22
	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,22
	   do 52 r=1,22
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,22
32	    pv(r)=tvn(uvs1,uvs2,uvs3,uvr(1,r),uvr(2,r),uvr(3,r))
	   pmn=10000.
           nrm=1
	   do 42 r=1,22
	      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,22
	   do 53 r=1,22
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,22
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,22
	      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(969),pr(969),484)
	rbi=sqrt(chisq/484.)
	call cchi(ob,pr,1452)

	RETURN
	entry pp(is,rmsd)
	write(7,117)is,rmsd,rbi,uar,uvr,uas,uvs
117	format('Subject 'i2' rmsd,rbi= '2f7.4,22(/' uar = '3f7.4),22(/' uvr = '3f7.4)
&                           ,22(/' uas = '3f7.4),22(/' 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

	subroutine rrn
	common /rnm/ mp,rr(1000000)
	mp=1
	return
	end

	subroutine crn
	common /rnm/ mp,rr(1000000)
	call srand(12345)
	do 10 j=1,1000000
10      rr(j)=crnorm(kk)
	mp=1
	return
	end

	function rnorm(n)
	common /rnm/ mp,rr(1000000)
	rnorm=rr(mp)
	mp=mp+1
	return
	end

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