(***********************************************************************
This file was generated automatically by the Mathematica front end.
It contains Initialization cells from a Notebook file, which typically
will have the same name as this file except ending in ".nb" instead of
".m".

This file is intended to be loaded into the Mathematica kernel using
the package loading commands Get or Needs.  Doing so is equivalent to
using the Evaluate Initialiation Cells menu command in the front end.

DO NOT EDIT THIS FILE.  This entire file is regenerated automatically 
each time the parent Notebook file is saved in the Mathematica front end.
Any changes you make to this file will be overwritten.
***********************************************************************)


(* Tilte: 
    Nonlinear Regression Analysis for Segregation Distortion in F2 Intercross 
      by OneOneZOne model *)

(* :Context: Application`OneOneZOne`   *)

(* :Name:OneOneZOne` *)

(* :Author: Yoshiaki Harushima *)

BeginPackage["SegAnalysis`OneOneZOne`", "SegAnalysis`jaco11Z11`", 
  "SegAnalysis`jaco11Z12`", "SegAnalysis`jaco11Z13`", 
  "SegAnalysis`jaco11Z14`",
	(* needed for graph *)
	"Graphics`MultipleListPlot`" ]

(*    *)
OneOneZOne::usage = 
  "OneOneZOne[dataA, dataB, dataH, renge, guess, {opts}] searches for \
positions and strengths of the segregation distortion genes by a \
least-squares fitting to lists of data within the renge according to the \
modele.  The model assumes two gametephyte genes and a zygotic viability gene \
to explain the marker segregation distortions on the chromosome.  The guess \
should be {{gametephyte gene position1, transmittance of A gamete}, \
{gametephyte gene position2, transmittance of A gamete}, {gene position, \
relative viability of B homozygote to A homozygote, relative viability of H \
homozygote to A homozygote}}.  The dataA, dataB and dataH can have the form \
{{x1, A1}, {x2, A2}, ...},  {{x1, B1}, {x2, B2}, ...},  and {{x1, H1}, {x2, \
H2}, ...}.  "

Options[OneOneZOne] = 
  Sort[Join[
		{AccuracyGoal -> Automatic, 
     Gradient -> Automatic, 
		MaxIterations -> 30, 
     PrecisionGoal -> Automatic, 
		Weights -> Equal, 
     WorkingPrecision -> $MachinePrecision,
		FixedParameters ->{}}, 
    Options[LinearSolve]]]

OneOneZOneA::usage =
  "OneOneZOneA[x, x1, t1, x2, t2, xv1, vb1, vh1] gives the expected fraction \
of A genotype at x, when gametophyte genes at x1 and x2 transmits A genotype \
by t1 and t2 on the different gamete and the zygotic vaibility gene at xv1 \
affects each viability of zygote.  The relateve viabilities of B homozygote \
and heterozygote to the A homozygote at xv1 are vb1 and vh1, respectively."

OneOneZOneB::usage =
  "OneOneZOneB[x, x1, t1, x2, t2, xv1, vb1, vh1] gives the expected fraction \
of B genotype at x, when gametophyte genes at x1 and x2 transmits A genotype \
by t1 and t2 on the different gamete and the zygotic vaibility gene at xv1 \
affects each viability of zygote.  The relateve viabilities of B homozygote \
and heterozygote to the A homozygote at xv1 are vb1 and vh1, respectively."

OneOneZOneH::usage =
  "OneOneZOneH[x, x1, t1, x2, t2, xv1, vb1, vh1] gives the expected fraction \
of H genotype at x, when gametophyte genes at x1 and x2 transmits A genotype \
by t1 and t2 on the different gamete and the zygotic vaibility gene at xv1 \
affects each viability of zygote.  The relateve viabilities of B homozygote \
and heterozygote to the A homozygote at xv1 are vb1 and vh1, respectively."

Begin["`Private`"]

(* === Warning messages ====  *)

OneOneZOne::baddataA = 
  "The dataA argument of OneOneZOne must be a vector."

OneOneZOne::baddataB = 
  "The dataB argument of OneOneZOne must be a vector."

OneOneZOne::baddataH = 
  "The dataH argument of OneOneZOne must be a vector."

OneOneZOne::badrenge = 
  "The renge argument of OneOneZOne must be `1` < `2`."

OneOneZOne::degrees = 
  "Warning: The data set is smaller than the number of parameters."

OneOneZOne::lmincrease = 
  "Warning: after `` decreases in the step size, NonlinearFit failed to find \
a new parameter estimate that will decrease the sum of squares."

OneOneZOne::badits =
"Warning: MaxIterations -> `1` is not set to a positive integer;
setting to 30."

OneOneZOne::badwork =
"Warning: The value `1` given for the WorkingPrecision is not an integer or \
is
an integer less than $MachinePrecision.  Setting to $MachinePrecision."

OneOneZOne::badprec =
"Warning: The value `1` given for the PrecisionGoal is not a positive integer
or Automatic.  Setting to $MachinePrecision - 10."

OneOneZOne::badacc =
"Warning: The value `1` given for the AccuracyGoal is not a positive integer
or Automatic.  Setting to $MachinePrecision - 10."

(* ===  OneOneZOne fatal error messages === *)

OneOneZOne::fitfail = "The fitting algorithm failed."

OneOneZOne::badguess = 
  "The guess should be {{gametephyte gene position1, transmittance of A \
gamete}, {gametephyte gene position2, transmittance of A gamete}, {gene \
position, relative viability of B homozygote to A homozygote, relative \
viability of H homozygote to A homozygote}}. Your guess is ``"
	
(* =========  OneOneZOne ===========  *)

(* OneOneZOne arguments;
	   dataA - matrix of numbers;
	   dataB - matrix of numbers;
	   dataH - matrix of numbers;
	   renge - the renge to fit modele equations;
	   guess - the initial guess;
	   opts - the options
		*)

OneOneZOne[dataA_List, dataB_List, dataH_List, renge_List, 
    guess_List, (opts___)?OptionQ]:=Module[{result, check}, 
    (
			result 
		)			/; (check = Check[
					 Module[{irenge, iguess}, 
					irenge = N[renge]; iguess = N[guess];
          result = oneonezone[dataA, dataB, dataH,irenge, iguess, opts]; 
           result =!= $Failed
						], False,OneOneZOne::fitfail];
				 check)
	] /; If[MatrixQ[N[dataA], NumberQ] || VectorQ[N[dataA], NumberQ], True, 
     Message[OneOneZOne::baddataA]; False] && 
    If[MatrixQ[N[dataB], NumberQ] || VectorQ[N[dataB], NumberQ], True, 
     Message[OneOneZOne::baddataB]; False] && 
    If[MatrixQ[N[dataH], NumberQ] || VectorQ[N[dataH], NumberQ], True, 
     Message[OneOneZOne::baddataH]; False]

oneonezone[dataA_, dataB_, dataH_, renge_, guess_, opts___] := 
  Module[{fixparams, maxits, workprec, accgoal, precgoal, modulus, zerotest, 
      minrenge, maxrenge, rdataA, rdataB, rdataH, ptsA, ptsB, ptsH, 
      responseA, responseB, responseH, prec, nptsA, nptsB, nptsH, x1, t1, x2, 
      t2, xv1, vb1, vh1, nresponseA, nresponseB, nresponseH, nx1, nt1, nx2, 
      nt2, nxv1, nvb1, nvh1},
		(* Get options *) 
   {fixparams, maxits, workprec, accgoal, precgoal, modulus, zerotest} = 
     {FixedParameters, MaxIterations, WorkingPrecision, AccuracyGoal, 
            PrecisionGoal, Modulus, ZeroTest} /. 
				{opts} /.  Options[OneOneZOne]; 
		(* Get and check renge *)
		 {minrenge, maxrenge} = renge;
    If[!(minrenge ===-Infinity || NumberQ[minrenge]) || 
      !(maxrenge === Infinity || NumberQ[maxrenge]) || 
      minrenge >= maxrenge, 
     Message[OneOneZOne::badrenge, minrenge, maxrenge]; Return[$Failed]];
		(*  - check guess *)
		{{x1, t1}, {x2, t2}, {xv1, vb1, vh1}} = guess;
		If[MatrixQ[guess, !NumberQ], Message[OneOneZOne::badguess, guess]; 
      Return[$Failed]];
		(* Select data by renge *)
		rdataA = Select[dataA, minrenge < #1[[1]] < maxrenge & ];
		rdataB = Select[dataB, minrenge < #1[[1]] < maxrenge & ];
		rdataH = Select[dataH, minrenge < #1[[1]] < maxrenge & ];
		(*  - assume univariate response *)
		responseA = Last /@ rdataA;
		ptsA = Flatten[(Take[#1, 1] & ) /@ rdataA];
		responseB = Last /@ rdataB;
		ptsB = Flatten[(Take[#1, 1] & ) /@ rdataB];
		responseH = Last /@ rdataH;
		ptsH =  Flatten[(Take[#1, 1] & ) /@ rdataH];
		(* - check degrees of freedom *)
		If[ Length[ptsA] + Length[ptsB] + Length[ptsH] < 7-Length[fixparams], 
      Message[OneOneZOne::degrees]];
		(*    *)
		If[ !(IntegerQ[maxits] && Positive[maxits]), 
     Message[OneOneZOne::badits, maxits]; maxits = 30];
		If[ !(IntegerQ[workprec] && 
        workprec >= $MachinePrecision), 
     Message[OneOneZOne::badwork, workprec]; 
      workprec = $MachinePrecision];
		If[precgoal === Automatic, precgoal = workprec - 10];
		If[accgoal === Automatic, accgoal = workprec - 10];
		If[ !(IntegerQ[precgoal] && Positive[precgoal]), 
     Message[OneOneZOne::badprec, precgoal]; 
      precgoal = $MachinePrecision - 10];
		If[ !(IntegerQ[accgoal] && Positive[accgoal]), 
     Message[OneOneZOne::badacc, accgoal]; 
      accgoal = $MachinePrecision - 10];
		(* - set working precision of inputs *)
		(* Note: 
        need the precision of the argument function to be at least as large 
          as WorkingPrecision *)
		nptsA = N[ptsA, workprec];
		nresponseA = N[responseA, workprec];
		nptsB = N[ptsB, workprec];
		nresponseB = N[responseB, workprec];
		nptsH = N[ptsH, workprec];
		nresponseH = N[responseH, workprec];
		(* - call proper fitting routine *)
		nx1 = N[x1, workprec]; nt1 = N[t1, workprec];
		nx2 = N[x2, workprec]; nt2 = N[t2, workprec];
		nxv1 = N[xv1, workprec];
		nvb1 = N[vb1, workprec]; nvh1 = N[vh1, workprec];
		out = mnlr[nptsA, nptsB, nptsH, nresponseA, nresponseB, nresponseH, nx1, 
        nt1, nx2, nt2, nxv1, nvb1, nvh1, fixparams, maxits, workprec, 
        precgoal, accgoal, modulus, zerotest];
		If[out =!= $Failed, Return[MakeOutputList[out]]];
]

(* =========    OneOneZOne     =========== *)
(* This uses Guess-Newton method for finding the least-
    squares solution for two gametophyte genes and a zygotic viability gene 
      model  *)

mnlr[nptsA_, nptsB_, nptsH_, nresponseA_, nresponseB_, nresponseH_, x1_, t1_, 
    x2_, t2_, xv1_, vb1_, vh1_, fp_, maxits_, wp_, pg_, ag_, mdls_, zt_] := 
  Module[{predA, predB, predH, chilistA, chilistB, chilistH, tmpchilistA, 
      tmpchilistB, tmpchilistH, tmpchilist, tmpnptsA, tmpnptsB, tmpnptsH, 
      chisq, derx1, dert1, derx2, dert2, derxv1, dervb1, dervh1, der, alpha, 
      beta, chisq0, newchi, tmpchi, its, accDelta, NAcc, precDelta, NPrec, 
      invaliddigits, accflag0, accflag1 = False, accflag2 = False, precflag0, 
      precflag1 = False, precflag2 = False, ovflag = False, guessx1 = x1, 
      guesst1 = t1, guessx2 = x2, guesst2 = t2, guessxv1 = xv1,
      guessvb1 = vb1, guessvh1 = vh1, oldx1, oldt1, oldx2, oldt2, oldxv1, 
      oldvb1, oldvh1, delta, deltax1, deltat1, deltax2, deltat2, deltaxv1, 
      deltavb1, deltavh1, tmpx1, tmpt1, tmpx2, tmpt2, tmpxv1, tmpvb1, tmpvh1, 
      dits, dmaxits, newpredA, newpredB, newpredH, newchilistA, newchilistB, 
      newchilistH, accgoal = N[10^(-ag)], precgoal = N[10^(-pg)], paramerr, 
      errorDOF},
		
		(* - calculation chi-squars *)
		{{chilistA, chilistB, chilistH}, chisq} = 
      calcchi[{nptsA, nptsB, nptsH}, {nresponseA, nresponseB,  nresponseH}, {
          guessx1, guesst1, guessx2, guesst2, guessxv1,guessvb1, guessvh1}, 
        wp];
		tmpchi = chisq; newchi = -1; its = 1;
		oldx1 = guessx1; oldt1 = guesst1;
		oldx2 = guessx2; oldt2 = guesst2;  
		oldxv1 = guessxv1;
		oldvb1 = guessvb1; oldvh1 = guessvh1;
		

		(* main loop *)
		ovflag = Check[
				While[(
						(accDelta = Abs[chisq - newchi];
							NAcc = -Floor[Log[10, accDelta]] - 1;
							precDelta = Abs[chisq - newchi]/(newchi + chisq)/2;
							NPrec = -Floor[Log[10, precDelta]] - 1;
							invaliddigits = 
                  TrueQ[NAcc =!= Infinity && Accuracy[accDelta] < NAcc && 
                      NPrec =!= Infinity && Precision[precDelta] < Nprec];
							!invaliddigits ) &&
							( accflag0 = ( 0 <= accDelta < accgoal) ;
								precflag0 =( 0 <= precDelta < precgoal);
								!((accflag0 && accflag1 && accflag2) || (
                        precflag0 && precflag1 && precflag2))  ) && 
              its <= maxits),
					
					{accflag2, precflag2} = {accflag1, precflag1};
					{accflag1, precflag1} = {accflag0, precflag0};
					chisq = tmpchi;
					oldx1 = guessx1; oldt1 = guesst1;
					oldx2 = guessx2; oldt2 = guesst2;
					oldxv1 = guessxv1; 
					oldvb1 = guessvb1; oldvh1 = guessvh1;
					
					(* === find deltax1, deltat1, deltax2, deltat2, deltaxv1, 
            deltavb1 and deltavh1 by LinearSolve normalmatrix ==== *)
					
					(* === delete not differentiable position ==== *)
					tmpchilistA = Delete[chilistA, Position[nptsA, guessx1]];
					tmpnptsA = Delete[nptsA, Position[nptsA, guessx1]];
					tmpchilistB =Delete[chilistB, Position[nptsB, guessx1]];
					tmpnptsB =Delete[nptsB, Position[nptsB, guessx1]];
					tmpchilistH =Delete[chilistH, Position[nptsH, guessx1]];
					tmpnptsH =Delete[nptsH, Position[nptsH, guessx1]];
					tmpchilistA = Delete[tmpchilistA, Position[tmpnptsA, guessx2]];
					tmpnptsA = Delete[tmpnptsA, Position[tmpnptsA, guessx2]];
					tmpchilistB =Delete[tmpchilistB, Position[tmpnptsB, guessx2]];
					tmpnptsB =Delete[tmpnptsB, Position[tmpnptsB, guessx2]];
					tmpchilistH =Delete[tmpchilistH, Position[tmpnptsH, guessx2]];
					tmpnptsH =Delete[tmpnptsH, Position[tmpnptsH, guessx2]]; 
					tmpchilistA = Delete[tmpchilistA, Position[tmpnptsA, guessxv1]];
					tmpnptsA =Delete[tmpnptsA, Position[tmpnptsA, guessxv1]];
					tmpchilistB = Delete[tmpchilistB, Position[tmpnptsB, guessxv1]];
					tmpnptsB =Delete[tmpnptsB, Position[tmpnptsB, guessxv1]];
					tmpchilistH = Delete[tmpchilistH, Position[tmpnptsH, guessxv1]];
					tmpnptsH =Delete[tmpnptsH, Position[tmpnptsH, guessxv1]];
					tmpchilist = Join[tmpchilistA , tmpchilistB, tmpchilistH];
					(* calculate jacobian *)
					
					der={};
					delta={};
					If[!(MemberQ[fp, "x1"]),
            derx1 =SegAnalysis`jaco11Z11`Private`a[guessx1, guesst1, guessx2, 
                guesst2, guessxv1, guessvb1, guessvh1, tmpnptsA, tmpnptsB, 
                tmpnptsH, wp];
						der=Append[der,derx1]];
					If[!(MemberQ[fp, "t1"]),
						dert1 =
              SegAnalysis`jaco11Z11`Private`b[guessx1, guesst1, guessx2, 
                guesst2, guessxv1, guessvb1, guessvh1, tmpnptsA, tmpnptsB, 
                tmpnptsH, wp];
						der=Append[der,dert1]];
					If[!(MemberQ[fp, "x2"]),
            derx2 =SegAnalysis`jaco11Z12`Private`a[guessx1, guesst1, guessx2, 
                guesst2, guessxv1, guessvb1, guessvh1, tmpnptsA, tmpnptsB, 
                tmpnptsH, wp];
						der=Append[der,derx2]];
					If[!(MemberQ[fp, "t2"]),
						dert2 =
              SegAnalysis`jaco11Z12`Private`b[guessx1, guesst1, guessx2, 
                guesst2, guessxv1, guessvb1, guessvh1, tmpnptsA, tmpnptsB, 
                tmpnptsH, wp];
						der=Append[der,dert2]];
					If[!(MemberQ[fp, "xv1"]),
						derxv1 =
              SegAnalysis`jaco11Z13`Private`a[guessx1, guesst1, guessx2, 
                guesst2, guessxv1, guessvb1, guessvh1, tmpnptsA, tmpnptsB, 
                tmpnptsH, wp];
						der=Append[der,derxv1]];
					If[!(MemberQ[fp, "vb1"]),
						dervb1 =
              SegAnalysis`jaco11Z13`Private`b[guessx1, guesst1, guessx2, 
                guesst2, guessxv1, guessvb1, guessvh1, tmpnptsA, tmpnptsB, 
                tmpnptsH, wp];
						der=Append[der,dervb1]];
					If[!(MemberQ[fp, "vh1"]),
						dervh1 =
              SegAnalysis`jaco11Z14`Private`a[guessx1, guesst1, guessx2, 
                guesst2, guessxv1, guessvb1, guessvh1, tmpnptsA, tmpnptsB, 
                tmpnptsH, wp];
						der=Append[der,dervh1]];
					
					(* normal equation: alpha . deltaguess === beta   *)
					{alpha, beta} = 
            N[{der . Transpose[der], Map[tmpchilist . #1 &, der]}, wp];
					
					delta = N[LinearSolve[alpha, beta], wp];
					If[!(MemberQ[fp, "x1"]),deltax1=First[delta];delta=Rest[delta], 
            deltax1=0];
					If[!(MemberQ[fp, "t1"]),deltat1=First[delta]; delta=Rest[delta],  
            deltat1=0];
					If[!(MemberQ[fp, "x2"]),deltax2=First[delta];delta=Rest[delta], 
            deltax2=0];
					If[!(MemberQ[fp, "t2"]),deltat2=First[delta];delta=Rest[delta],  
            deltat2=0];
					If[!(MemberQ[fp, "xv1"]),deltaxv1=First[delta];delta=Rest[delta], 
            deltaxv1=0];
					If[!(MemberQ[fp, "vb1"]),deltavb1=First[delta]; delta=Rest[delta],  
            deltavb1=0];
					If[!(MemberQ[fp, "vh1"]),deltavh1=First[delta],  deltavh1=0];

					
					{tmpx1, tmpt1, tmpx2, tmpt2, tmpxv1,tmpvb1, tmpvh1} = 
            N[{oldx1, oldt1, oldx2, oldt2, oldxv1, oldvb1, oldvh1} +  {
                  deltax1, deltat1, deltax2, deltat2, deltaxv1, deltavb1, 
                  deltavh1}, wp];
					
					{{newchilistA, newchilistB, newchilistH}, newchi} = 
            calcchi[{nptsA, nptsB, nptsH}, {nresponseA, nresponseB, 
                nresponseH}, {tmpx1, tmpt1, tmpx2, tmpt2, tmpxv1,tmpvb1, 
                tmpvh1}, wp];
					
					(* find newguess for smaller chisq  *)
					dits = 1; dmaxits = 15;
					While[newchi - chisq >= 0 && dits <= dmaxits,
						{tmpx1, tmpt1, tmpx2, tmpt2, tmpxv1,tmpvb1, tmpvh1} = 
              N[{oldx1, oldt1, oldx2, oldt2, oldxv1, oldvb1, 
                    oldvh1} + (2^(-1*dits))*{deltax1, deltat1, deltax2, 
                      deltat2, deltaxv1, deltavb1, deltavh1}, wp];
						{{newchilistA, newchilistB, newchilistH}, newchi} = 
              calcchi[{nptsA, nptsB, nptsH}, {nresponseA, nresponseB, 
                  nresponseH}, {tmpx1, tmpt1, tmpx2, tmpt2, tmpxv1,tmpvb1, 
                  tmpvh1}, wp];
						dits++];
					If[dits > dmaxits, Message[OneOneZOne::lmincrease, dmaxits]];
					tmpchi = newchi;
					{guessx1, guesst1, guessx2, guesst2, guessxv1, guessvb1, guessvh1} = {
              tmpx1, tmpt1, tmpx2, tmpt2, tmpxv1,tmpvb1, tmpvh1};
					chilistA = newchilistA;
					chilistB = newchilistB;
					chilistH = newchilistH;
					printprogress[its, 
            newchi, {guessx1, guesst1, guessx2, guesst2, guessxv1, guessvb1, 
              guessvh1}];
					its++], (* end of While *)
				True, General::ovfl];  (* end of Check *)
		(* Adjust the accuracy and precision of chisq *)
		If[Accuracy[newchi] > NAcc, newchi = SetAccuracy[newchi, NAcc]];
		If[Precision[newchi] > NPrec, newchi = SetPrecision[newchi, NPrec]];
		errorDOF = 
      Length[chilistA] + Length[chilistB] + Length[chilistH] - 7+Length[fp];
		paramerr = Sqrt[ ( newchi/errorDOF) DiagonalElements[Inverse[alpha]]];
		{newchi, {guessx1, guesst1, guessx2, guesst2, guessxv1, guessvb1, 
        guessvh1}, paramerr, errorDOF, fp}
	]  (* end OneOneZOne *)

(*   *)
OneOneZOneA[x_, x1_, t1_, x2_, t2_, xv1_, vb1_, vh1_] := Which[
		x<= x1 && x <= x2 && x <= xv1,
		(8*Cosh[2*(x - xv1)]^2*(
            2*t1*Cosh[2*(x - x1)] + (-1 + 2*t1)*Sinh[2*(x - x1)])*(
            2*t2*Cosh[2*(x - x2)] + (-1 + 2*t2)*Sinh[2*(x - x2)])*(
            1 + ((1 + vb1 - 2*vh1)*Tanh[2*(x - xv1)]^2)/4 + (-1 + vh1)*
                Tanh[2*(-x + xv1)]))/(-((
                  Cosh[2*(x1 - x2)] + (-1 + 4*t1 + 4*t2 - 8*t1*t2)*
                      Cosh[4*x - 2*(x1 + x2)] - 
                    2*(-1 + 2*t1)*(-1 + 2*t2)*Sinh[4*x - 2*(x1 + x2)])*(-1 - 
                    vb1 - 2*vh1 + (1 + vb1 - 6*vh1)*Cosh[4*(x - xv1)] + 
                    2*(1 + vb1 - 2*vh1)*Sinh[4*(x - xv1)])) + (
              2*(-1 + t1)*Cosh[2*(x - x1)] + (-1 + 2*t1)*Sinh[2*(x - x1)])*(
              2*(-1 + t2)*Cosh[2*(x - x2)] + (-1 + 2*t2)*
                  Sinh[2*(x - x2)])*(-1 + 3*vb1 + 
                2*vh1 + (1 + 5*vb1 - 2*vh1)*Cosh[4*(x - xv1)] + 
                4*(vb1 - vh1)*Sinh[4*(x - xv1)]) + (
              2*t1*Cosh[2*(x - x1)] + (-1 + 2*t1)*Sinh[2*(x - x1)])*(
              2*t2*Cosh[2*(x - x2)] + (-1 + 2*t2)*Sinh[2*(x - x2)])*(
              3 - vb1 + 2*vh1 + (5 + vb1 - 2*vh1)*Cosh[4*(x - xv1)] - 
                4*(-1 + vh1)*Sinh[4*(x - xv1)])),
		x<= x1 && x <= x2 && xv1 < x,
		(8*Cosh[2*(x - xv1)]^2*(
            2*t1*Cosh[2*(x - x1)] + (-1 + 2*t1)*Sinh[2*(x - x1)])*(
            2*t2*Cosh[2*(x - x2)] + (-1 + 2*t2)*Sinh[2*(x - x2)])*(
            1 + (-1 + vh1)*
                Tanh[2*x - 2*xv1] + ((1 + vb1 - 2*vh1)*Tanh[2*x - 2*xv1]^2)/
                4))/(-((
                  Cosh[2*(x1 - x2)] + (-1 + 4*t1 + 4*t2 - 8*t1*t2)*
                      Cosh[4*x - 2*(x1 + x2)] - 
                    2*(-1 + 2*t1)*(-1 + 2*t2)*Sinh[4*x - 2*(x1 + x2)])*(-1 - 
                    vb1 - 2*vh1 + (1 + vb1 - 6*vh1)*Cosh[4*(x - xv1)] - 
                    2*(1 + vb1 - 2*vh1)*Sinh[4*(x - xv1)])) + (
              2*(-1 + t1)*Cosh[2*(x - x1)] + (-1 + 2*t1)*Sinh[2*(x - x1)])*(
              2*(-1 + t2)*Cosh[2*(x - x2)] + (-1 + 2*t2)*
                  Sinh[2*(x - x2)])*(-1 + 3*vb1 + 
                2*vh1 + (1 + 5*vb1 - 2*vh1)*Cosh[4*(x - xv1)] - 
                4*(vb1 - vh1)*Sinh[4*(x - xv1)]) + (
              2*t1*Cosh[2*(x - x1)] + (-1 + 2*t1)*Sinh[2*(x - x1)])*(
              2*t2*Cosh[2*(x - x2)] + (-1 + 2*t2)*Sinh[2*(x - x2)])*(
              3 - vb1 + 2*vh1 + (5 + vb1 - 2*vh1)*Cosh[4*(x - xv1)] + 
                4*(-1 + vh1)*Sinh[4*(x - xv1)])),
		x<= x1 && x2 < x && x <= xv1,
		(8*Cosh[2*(x - xv1)]^2*(
            2*t1*Cosh[2*(x - x1)] + (-1 + 2*t1)*Sinh[2*(x - x1)])*(
            2*t2*Cosh[2*(x - x2)] + (1 - 2*t2)*Sinh[2*(x - x2)])*(
            1 + ((1 + vb1 - 2*vh1)*Tanh[2*(x - xv1)]^2)/4 + (-1 + vh1)*
                Tanh[2*(-x + xv1)]))/(-(((-1 + t1*(4 - 8*t2) + 4*t2)*
                      Cosh[2*(x1 - x2)] + Cosh[4*x - 2*(x1 + x2)] + 
                    2*(-1 + 2*t1)*(-1 + 2*t2)*Sinh[2*(x1 - x2)])*(-1 - vb1 - 
                    2*vh1 + (1 + vb1 - 6*vh1)*Cosh[4*(x - xv1)] + 
                    2*(1 + vb1 - 2*vh1)*Sinh[4*(x - xv1)])) + (
              2*(-1 + t1)*Cosh[2*(x - x1)] + (-1 + 2*t1)*Sinh[2*(x - x1)])*(
              2*(-1 + t2)*Cosh[2*(x - x2)] + (1 - 2*t2)*
                  Sinh[2*(x - x2)])*(-1 + 3*vb1 + 
                2*vh1 + (1 + 5*vb1 - 2*vh1)*Cosh[4*(x - xv1)] + 
                4*(vb1 - vh1)*Sinh[4*(x - xv1)]) + (
              2*t1*Cosh[2*(x - x1)] + (-1 + 2*t1)*Sinh[2*(x - x1)])*(
              2*t2*Cosh[2*(x - x2)] + (1 - 2*t2)*Sinh[2*(x - x2)])*(
              3 - vb1 + 2*vh1 + (5 + vb1 - 2*vh1)*Cosh[4*(x - xv1)] - 
                4*(-1 + vh1)*Sinh[4*(x - xv1)])),
		x<= x1 && x2 < x && xv1 < x,
		(8*Cosh[2*(x - xv1)]^2*(
            2*t1*Cosh[2*(x - x1)] + (-1 + 2*t1)*Sinh[2*(x - x1)])*(
            2*t2*Cosh[2*(x - x2)] + (1 - 2*t2)*Sinh[2*(x - x2)])*(
            1 + (-1 + vh1)*
                Tanh[2*x - 2*xv1] + ((1 + vb1 - 2*vh1)*Tanh[2*x - 2*xv1]^2)/
                4))/(-(((-1 + t1*(4 - 8*t2) + 4*t2)*Cosh[2*(x1 - x2)] + 
                    Cosh[4*x - 2*(x1 + x2)] + 
                    2*(-1 + 2*t1)*(-1 + 2*t2)*Sinh[2*(x1 - x2)])*(-1 - vb1 - 
                    2*vh1 + (1 + vb1 - 6*vh1)*Cosh[4*(x - xv1)] - 
                    2*(1 + vb1 - 2*vh1)*Sinh[4*(x - xv1)])) + (
              2*(-1 + t1)*Cosh[2*(x - x1)] + (-1 + 2*t1)*Sinh[2*(x - x1)])*(
              2*(-1 + t2)*Cosh[2*(x - x2)] + (1 - 2*t2)*
                  Sinh[2*(x - x2)])*(-1 + 3*vb1 + 
                2*vh1 + (1 + 5*vb1 - 2*vh1)*Cosh[4*(x - xv1)] - 
                4*(vb1 - vh1)*Sinh[4*(x - xv1)]) + (
              2*t1*Cosh[2*(x - x1)] + (-1 + 2*t1)*Sinh[2*(x - x1)])*(
              2*t2*Cosh[2*(x - x2)] + (1 - 2*t2)*Sinh[2*(x - x2)])*(
              3 - vb1 + 2*vh1 + (5 + vb1 - 2*vh1)*Cosh[4*(x - xv1)] + 
                4*(-1 + vh1)*Sinh[4*(x - xv1)])),
		x1< x && x <= x2 && x <= xv1,
		(8*Cosh[2*(x - xv1)]^2*(
            2*t1*Cosh[2*(x - x1)] + (1 - 2*t1)*Sinh[2*(x - x1)])*(
            2*t2*Cosh[2*(x - x2)] + (-1 + 2*t2)*Sinh[2*(x - x2)])*(
            1 + ((1 + vb1 - 2*vh1)*Tanh[2*(x - xv1)]^2)/4 + (-1 + vh1)*
                Tanh[2*(-x + xv1)]))/(-(((-1 + t1*(4 - 8*t2) + 4*t2)*
                      Cosh[2*(x1 - x2)] + Cosh[4*x - 2*(x1 + x2)] - 
                    2*(-1 + 2*t1)*(-1 + 2*t2)*Sinh[2*(x1 - x2)])*(-1 - vb1 - 
                    2*vh1 + (1 + vb1 - 6*vh1)*Cosh[4*(x - xv1)] + 
                    2*(1 + vb1 - 2*vh1)*Sinh[4*(x - xv1)])) + (
              2*(-1 + t1)*Cosh[2*(x - x1)] + (1 - 2*t1)*Sinh[2*(x - x1)])*(
              2*(-1 + t2)*Cosh[2*(x - x2)] + (-1 + 2*t2)*
                  Sinh[2*(x - x2)])*(-1 + 3*vb1 + 
                2*vh1 + (1 + 5*vb1 - 2*vh1)*Cosh[4*(x - xv1)] + 
                4*(vb1 - vh1)*Sinh[4*(x - xv1)]) + (
              2*t1*Cosh[2*(x - x1)] + (1 - 2*t1)*Sinh[2*(x - x1)])*(
              2*t2*Cosh[2*(x - x2)] + (-1 + 2*t2)*Sinh[2*(x - x2)])*(
              3 - vb1 + 2*vh1 + (5 + vb1 - 2*vh1)*Cosh[4*(x - xv1)] - 
                4*(-1 + vh1)*Sinh[4*(x - xv1)])),
		x1< x && x <= x2 && xv1 < x,
		(8*Cosh[2*(x - xv1)]^2*(
            2*t1*Cosh[2*(x - x1)] + (1 - 2*t1)*Sinh[2*(x - x1)])*(
            2*t2*Cosh[2*(x - x2)] + (-1 + 2*t2)*Sinh[2*(x - x2)])*(
            1 + (-1 + vh1)*
                Tanh[2*x - 2*xv1] + ((1 + vb1 - 2*vh1)*Tanh[2*x - 2*xv1]^2)/
                4))/(-(((-1 + t1*(4 - 8*t2) + 4*t2)*Cosh[2*(x1 - x2)] + 
                    Cosh[4*x - 2*(x1 + x2)] - 
                    2*(-1 + 2*t1)*(-1 + 2*t2)*Sinh[2*(x1 - x2)])*(-1 - vb1 - 
                    2*vh1 + (1 + vb1 - 6*vh1)*Cosh[4*(x - xv1)] - 
                    2*(1 + vb1 - 2*vh1)*Sinh[4*(x - xv1)])) + (
              2*(-1 + t1)*Cosh[2*(x - x1)] + (1 - 2*t1)*Sinh[2*(x - x1)])*(
              2*(-1 + t2)*Cosh[2*(x - x2)] + (-1 + 2*t2)*
                  Sinh[2*(x - x2)])*(-1 + 3*vb1 + 
                2*vh1 + (1 + 5*vb1 - 2*vh1)*Cosh[4*(x - xv1)] - 
                4*(vb1 - vh1)*Sinh[4*(x - xv1)]) + (
              2*t1*Cosh[2*(x - x1)] + (1 - 2*t1)*Sinh[2*(x - x1)])*(
              2*t2*Cosh[2*(x - x2)] + (-1 + 2*t2)*Sinh[2*(x - x2)])*(
              3 - vb1 + 2*vh1 + (5 + vb1 - 2*vh1)*Cosh[4*(x - xv1)] + 
                4*(-1 + vh1)*Sinh[4*(x - xv1)])),
		x1< x && x2 < x && x <= xv1,
		(8*Cosh[2*(x - xv1)]^2*(
            2*t1*Cosh[2*(x - x1)] + (1 - 2*t1)*Sinh[2*(x - x1)])*(
            2*t2*Cosh[2*(x - x2)] + (1 - 2*t2)*Sinh[2*(x - x2)])*(
            1 + ((1 + vb1 - 2*vh1)*Tanh[2*(x - xv1)]^2)/4 + (-1 + vh1)*
                Tanh[2*(-x + xv1)]))/(-((
                  Cosh[2*(x1 - x2)] + (-1 + 4*t1 + 4*t2 - 8*t1*t2)*
                      Cosh[4*x - 2*(x1 + x2)] + 
                    2*(-1 + 2*t1)*(-1 + 2*t2)*Sinh[4*x - 2*(x1 + x2)])*(-1 - 
                    vb1 - 2*vh1 + (1 + vb1 - 6*vh1)*Cosh[4*(x - xv1)] + 
                    2*(1 + vb1 - 2*vh1)*Sinh[4*(x - xv1)])) + (
              2*(-1 + t1)*Cosh[2*(x - x1)] + (1 - 2*t1)*Sinh[2*(x - x1)])*(
              2*(-1 + t2)*Cosh[2*(x - x2)] + (1 - 2*t2)*
                  Sinh[2*(x - x2)])*(-1 + 3*vb1 + 
                2*vh1 + (1 + 5*vb1 - 2*vh1)*Cosh[4*(x - xv1)] + 
                4*(vb1 - vh1)*Sinh[4*(x - xv1)]) + (
              2*t1*Cosh[2*(x - x1)] + (1 - 2*t1)*Sinh[2*(x - x1)])*(
              2*t2*Cosh[2*(x - x2)] + (1 - 2*t2)*Sinh[2*(x - x2)])*(
              3 - vb1 + 2*vh1 + (5 + vb1 - 2*vh1)*Cosh[4*(x - xv1)] - 
                4*(-1 + vh1)*Sinh[4*(x - xv1)])),
		x1< x && x2 < x && xv1 < x,
		(8*Cosh[2*(x - xv1)]^2*(
            2*t1*Cosh[2*(x - x1)] + (1 - 2*t1)*Sinh[2*(x - x1)])*(
            2*t2*Cosh[2*(x - x2)] + (1 - 2*t2)*Sinh[2*(x - x2)])*(
            1 + (-1 + vh1)*
                Tanh[2*x - 2*xv1] + ((1 + vb1 - 2*vh1)*Tanh[2*x - 2*xv1]^2)/
                4))/(-((
                  Cosh[2*(x1 - x2)] + (-1 + 4*t1 + 4*t2 - 8*t1*t2)*
                      Cosh[4*x - 2*(x1 + x2)] + 
                    2*(-1 + 2*t1)*(-1 + 2*t2)*Sinh[4*x - 2*(x1 + x2)])*(-1 - 
                    vb1 - 2*vh1 + (1 + vb1 - 6*vh1)*Cosh[4*(x - xv1)] - 
                    2*(1 + vb1 - 2*vh1)*Sinh[4*(x - xv1)])) + (
              2*(-1 + t1)*Cosh[2*(x - x1)] + (1 - 2*t1)*Sinh[2*(x - x1)])*(
              2*(-1 + t2)*Cosh[2*(x - x2)] + (1 - 2*t2)*
                  Sinh[2*(x - x2)])*(-1 + 3*vb1 + 
                2*vh1 + (1 + 5*vb1 - 2*vh1)*Cosh[4*(x - xv1)] - 
                4*(vb1 - vh1)*Sinh[4*(x - xv1)]) + (
              2*t1*Cosh[2*(x - x1)] + (1 - 2*t1)*Sinh[2*(x - x1)])*(
              2*t2*Cosh[2*(x - x2)] + (1 - 2*t2)*Sinh[2*(x - x2)])*(
              3 - vb1 + 2*vh1 + (5 + vb1 - 2*vh1)*Cosh[4*(x - xv1)] + 
                4*(-1 + vh1)*Sinh[4*(x - xv1)]))]
	
(*  *)
OneOneZOneB[x_, x1_, t1_, x2_, t2_, xv1_, vb1_, vh1_] := Which[
		x<= x1 && x <= x2 && x <= xv1,
		(2*Cosh[2*(x - xv1)]^2*(
            2*(-1 + t1)*Cosh[2*(x - x1)] + (-1 + 2*t1)*Sinh[2*(x - x1)])*(
            2*(-1 + t2)*Cosh[2*(x - x2)] + (-1 + 2*t2)*Sinh[2*(x - x2)])*(
            4*vb1*Tanh[2*(x - xv1)] + (1 + vb1 - 2*vh1)*Tanh[2*(x - xv1)]^2 + 
              4*(vb1 + vh1*Tanh[2*(-x + xv1)])))/(-((
                  Cosh[2*(x1 - x2)] + (-1 + 4*t1 + 4*t2 - 8*t1*t2)*
                      Cosh[4*x - 2*(x1 + x2)] - 
                    2*(-1 + 2*t1)*(-1 + 2*t2)*Sinh[4*x - 2*(x1 + x2)])*(-1 - 
                    vb1 - 2*vh1 + (1 + vb1 - 6*vh1)*Cosh[4*(x - xv1)] + 
                    2*(1 + vb1 - 2*vh1)*Sinh[4*(x - xv1)])) + (
              2*(-1 + t1)*Cosh[2*(x - x1)] + (-1 + 2*t1)*Sinh[2*(x - x1)])*(
              2*(-1 + t2)*Cosh[2*(x - x2)] + (-1 + 2*t2)*
                  Sinh[2*(x - x2)])*(-1 + 3*vb1 + 
                2*vh1 + (1 + 5*vb1 - 2*vh1)*Cosh[4*(x - xv1)] + 
                4*(vb1 - vh1)*Sinh[4*(x - xv1)]) + (
              2*t1*Cosh[2*(x - x1)] + (-1 + 2*t1)*Sinh[2*(x - x1)])*(
              2*t2*Cosh[2*(x - x2)] + (-1 + 2*t2)*Sinh[2*(x - x2)])*(
              3 - vb1 + 2*vh1 + (5 + vb1 - 2*vh1)*Cosh[4*(x - xv1)] - 
                4*(-1 + vh1)*Sinh[4*(x - xv1)])),
		x<= x1 && x <= x2 && xv1 < x,
		(2*Cosh[2*(x - xv1)]^2*(
            2*(-1 + t1)*Cosh[2*(x - x1)] + (-1 + 2*t1)*Sinh[2*(x - x1)])*(
            2*(-1 + t2)*Cosh[2*(x - x2)] + (-1 + 2*t2)*Sinh[2*(x - x2)])*(
            4*vb1 - 4*(vb1 - vh1)*Tanh[2*x - 2*xv1] + (1 + vb1 - 2*vh1)*
                Tanh[2*x - 2*xv1]^2))/(-((
                  Cosh[2*(x1 - x2)] + (-1 + 4*t1 + 4*t2 - 8*t1*t2)*
                      Cosh[4*x - 2*(x1 + x2)] - 
                    2*(-1 + 2*t1)*(-1 + 2*t2)*Sinh[4*x - 2*(x1 + x2)])*(-1 - 
                    vb1 - 2*vh1 + (1 + vb1 - 6*vh1)*Cosh[4*(x - xv1)] - 
                    2*(1 + vb1 - 2*vh1)*Sinh[4*(x - xv1)])) + (
              2*(-1 + t1)*Cosh[2*(x - x1)] + (-1 + 2*t1)*Sinh[2*(x - x1)])*(
              2*(-1 + t2)*Cosh[2*(x - x2)] + (-1 + 2*t2)*
                  Sinh[2*(x - x2)])*(-1 + 3*vb1 + 
                2*vh1 + (1 + 5*vb1 - 2*vh1)*Cosh[4*(x - xv1)] - 
                4*(vb1 - vh1)*Sinh[4*(x - xv1)]) + (
              2*t1*Cosh[2*(x - x1)] + (-1 + 2*t1)*Sinh[2*(x - x1)])*(
              2*t2*Cosh[2*(x - x2)] + (-1 + 2*t2)*Sinh[2*(x - x2)])*(
              3 - vb1 + 2*vh1 + (5 + vb1 - 2*vh1)*Cosh[4*(x - xv1)] + 
                4*(-1 + vh1)*Sinh[4*(x - xv1)])),
		x<= x1 && x2 < x && x <= xv1,
		(2*Cosh[2*(x - xv1)]^2*(
            2*(-1 + t1)*Cosh[2*(x - x1)] + (-1 + 2*t1)*Sinh[2*(x - x1)])*(
            2*(-1 + t2)*Cosh[2*(x - x2)] + (1 - 2*t2)*Sinh[2*(x - x2)])*(
            4*vb1*Tanh[2*(x - xv1)] + (1 + vb1 - 2*vh1)*Tanh[2*(x - xv1)]^2 + 
              4*(vb1 + 
                    vh1*Tanh[2*(-x + xv1)])))/(-(((-1 + t1*(4 - 8*t2) + 4*t2)*
                      Cosh[2*(x1 - x2)] + Cosh[4*x - 2*(x1 + x2)] + 
                    2*(-1 + 2*t1)*(-1 + 2*t2)*Sinh[2*(x1 - x2)])*(-1 - vb1 - 
                    2*vh1 + (1 + vb1 - 6*vh1)*Cosh[4*(x - xv1)] + 
                    2*(1 + vb1 - 2*vh1)*Sinh[4*(x - xv1)])) + (
              2*(-1 + t1)*Cosh[2*(x - x1)] + (-1 + 2*t1)*Sinh[2*(x - x1)])*(
              2*(-1 + t2)*Cosh[2*(x - x2)] + (1 - 2*t2)*
                  Sinh[2*(x - x2)])*(-1 + 3*vb1 + 
                2*vh1 + (1 + 5*vb1 - 2*vh1)*Cosh[4*(x - xv1)] + 
                4*(vb1 - vh1)*Sinh[4*(x - xv1)]) + (
              2*t1*Cosh[2*(x - x1)] + (-1 + 2*t1)*Sinh[2*(x - x1)])*(
              2*t2*Cosh[2*(x - x2)] + (1 - 2*t2)*Sinh[2*(x - x2)])*(
              3 - vb1 + 2*vh1 + (5 + vb1 - 2*vh1)*Cosh[4*(x - xv1)] - 
                4*(-1 + vh1)*Sinh[4*(x - xv1)])),
		x<= x1 && x2 < x && xv1 < x,
		(2*Cosh[2*(x - xv1)]^2*(
            2*(-1 + t1)*Cosh[2*(x - x1)] + (-1 + 2*t1)*Sinh[2*(x - x1)])*(
            2*(-1 + t2)*Cosh[2*(x - x2)] + (1 - 2*t2)*Sinh[2*(x - x2)])*(
            4*vb1 - 4*(vb1 - vh1)*Tanh[2*x - 2*xv1] + (1 + vb1 - 2*vh1)*
                Tanh[2*x - 2*xv1]^2))/(-(((-1 + t1*(4 - 8*t2) + 4*t2)*
                      Cosh[2*(x1 - x2)] + Cosh[4*x - 2*(x1 + x2)] + 
                    2*(-1 + 2*t1)*(-1 + 2*t2)*Sinh[2*(x1 - x2)])*(-1 - vb1 - 
                    2*vh1 + (1 + vb1 - 6*vh1)*Cosh[4*(x - xv1)] - 
                    2*(1 + vb1 - 2*vh1)*Sinh[4*(x - xv1)])) + (
              2*(-1 + t1)*Cosh[2*(x - x1)] + (-1 + 2*t1)*Sinh[2*(x - x1)])*(
              2*(-1 + t2)*Cosh[2*(x - x2)] + (1 - 2*t2)*
                  Sinh[2*(x - x2)])*(-1 + 3*vb1 + 
                2*vh1 + (1 + 5*vb1 - 2*vh1)*Cosh[4*(x - xv1)] - 
                4*(vb1 - vh1)*Sinh[4*(x - xv1)]) + (
              2*t1*Cosh[2*(x - x1)] + (-1 + 2*t1)*Sinh[2*(x - x1)])*(
              2*t2*Cosh[2*(x - x2)] + (1 - 2*t2)*Sinh[2*(x - x2)])*(
              3 - vb1 + 2*vh1 + (5 + vb1 - 2*vh1)*Cosh[4*(x - xv1)] + 
                4*(-1 + vh1)*Sinh[4*(x - xv1)])),
		x1< x && x <= x2 && x <= xv1,
		(2*Cosh[2*(x - xv1)]^2*(
            2*(-1 + t1)*Cosh[2*(x - x1)] + (1 - 2*t1)*Sinh[2*(x - x1)])*(
            2*(-1 + t2)*Cosh[2*(x - x2)] + (-1 + 2*t2)*Sinh[2*(x - x2)])*(
            4*vb1*Tanh[2*(x - xv1)] + (1 + vb1 - 2*vh1)*Tanh[2*(x - xv1)]^2 + 
              4*(vb1 + 
                    vh1*Tanh[2*(-x + xv1)])))/(-(((-1 + t1*(4 - 8*t2) + 4*t2)*
                      Cosh[2*(x1 - x2)] + Cosh[4*x - 2*(x1 + x2)] - 
                    2*(-1 + 2*t1)*(-1 + 2*t2)*Sinh[2*(x1 - x2)])*(-1 - vb1 - 
                    2*vh1 + (1 + vb1 - 6*vh1)*Cosh[4*(x - xv1)] + 
                    2*(1 + vb1 - 2*vh1)*Sinh[4*(x - xv1)])) + (
              2*(-1 + t1)*Cosh[2*(x - x1)] + (1 - 2*t1)*Sinh[2*(x - x1)])*(
              2*(-1 + t2)*Cosh[2*(x - x2)] + (-1 + 2*t2)*
                  Sinh[2*(x - x2)])*(-1 + 3*vb1 + 
                2*vh1 + (1 + 5*vb1 - 2*vh1)*Cosh[4*(x - xv1)] + 
                4*(vb1 - vh1)*Sinh[4*(x - xv1)]) + (
              2*t1*Cosh[2*(x - x1)] + (1 - 2*t1)*Sinh[2*(x - x1)])*(
              2*t2*Cosh[2*(x - x2)] + (-1 + 2*t2)*Sinh[2*(x - x2)])*(
              3 - vb1 + 2*vh1 + (5 + vb1 - 2*vh1)*Cosh[4*(x - xv1)] - 
                4*(-1 + vh1)*Sinh[4*(x - xv1)])),
		x1< x && x <= x2 && xv1 < x,
		(2*Cosh[2*(x - xv1)]^2*(
            2*(-1 + t1)*Cosh[2*(x - x1)] + (1 - 2*t1)*Sinh[2*(x - x1)])*(
            2*(-1 + t2)*Cosh[2*(x - x2)] + (-1 + 2*t2)*Sinh[2*(x - x2)])*(
            4*vb1 - 4*(vb1 - vh1)*Tanh[2*x - 2*xv1] + (1 + vb1 - 2*vh1)*
                Tanh[2*x - 2*xv1]^2))/(-(((-1 + t1*(4 - 8*t2) + 4*t2)*
                      Cosh[2*(x1 - x2)] + Cosh[4*x - 2*(x1 + x2)] - 
                    2*(-1 + 2*t1)*(-1 + 2*t2)*Sinh[2*(x1 - x2)])*(-1 - vb1 - 
                    2*vh1 + (1 + vb1 - 6*vh1)*Cosh[4*(x - xv1)] - 
                    2*(1 + vb1 - 2*vh1)*Sinh[4*(x - xv1)])) + (
              2*(-1 + t1)*Cosh[2*(x - x1)] + (1 - 2*t1)*Sinh[2*(x - x1)])*(
              2*(-1 + t2)*Cosh[2*(x - x2)] + (-1 + 2*t2)*
                  Sinh[2*(x - x2)])*(-1 + 3*vb1 + 
                2*vh1 + (1 + 5*vb1 - 2*vh1)*Cosh[4*(x - xv1)] - 
                4*(vb1 - vh1)*Sinh[4*(x - xv1)]) + (
              2*t1*Cosh[2*(x - x1)] + (1 - 2*t1)*Sinh[2*(x - x1)])*(
              2*t2*Cosh[2*(x - x2)] + (-1 + 2*t2)*Sinh[2*(x - x2)])*(
              3 - vb1 + 2*vh1 + (5 + vb1 - 2*vh1)*Cosh[4*(x - xv1)] + 
                4*(-1 + vh1)*Sinh[4*(x - xv1)])),
		x1< x && x2 < x && x <= xv1,
		(2*Cosh[2*(x - xv1)]^2*(
            2*(-1 + t1)*Cosh[2*(x - x1)] + (1 - 2*t1)*Sinh[2*(x - x1)])*(
            2*(-1 + t2)*Cosh[2*(x - x2)] + (1 - 2*t2)*Sinh[2*(x - x2)])*(
            4*vb1*Tanh[2*(x - xv1)] + (1 + vb1 - 2*vh1)*Tanh[2*(x - xv1)]^2 + 
              4*(vb1 + vh1*Tanh[2*(-x + xv1)])))/(-((
                  Cosh[2*(x1 - x2)] + (-1 + 4*t1 + 4*t2 - 8*t1*t2)*
                      Cosh[4*x - 2*(x1 + x2)] + 
                    2*(-1 + 2*t1)*(-1 + 2*t2)*Sinh[4*x - 2*(x1 + x2)])*(-1 - 
                    vb1 - 2*vh1 + (1 + vb1 - 6*vh1)*Cosh[4*(x - xv1)] + 
                    2*(1 + vb1 - 2*vh1)*Sinh[4*(x - xv1)])) + (
              2*(-1 + t1)*Cosh[2*(x - x1)] + (1 - 2*t1)*Sinh[2*(x - x1)])*(
              2*(-1 + t2)*Cosh[2*(x - x2)] + (1 - 2*t2)*
                  Sinh[2*(x - x2)])*(-1 + 3*vb1 + 
                2*vh1 + (1 + 5*vb1 - 2*vh1)*Cosh[4*(x - xv1)] + 
                4*(vb1 - vh1)*Sinh[4*(x - xv1)]) + (
              2*t1*Cosh[2*(x - x1)] + (1 - 2*t1)*Sinh[2*(x - x1)])*(
              2*t2*Cosh[2*(x - x2)] + (1 - 2*t2)*Sinh[2*(x - x2)])*(
              3 - vb1 + 2*vh1 + (5 + vb1 - 2*vh1)*Cosh[4*(x - xv1)] - 
                4*(-1 + vh1)*Sinh[4*(x - xv1)])),
		x1< x && x2 < x && xv1 < x,
		((2*(-1 + t1)*Cosh[2*(x - x1)] + (1 - 2*t1)*Sinh[2*(x - x1)])*(
            2*(-1 + t2)*Cosh[2*(x - x2)] + (1 - 2*t2)*Sinh[2*(x - x2)])*(-1 + 
              3*vb1 + 2*vh1 + (1 + 5*vb1 - 2*vh1)*Cosh[4*(x - xv1)] - 
              4*(vb1 - vh1)*
                Sinh[4*(x - xv1)]))/(-((
                  Cosh[2*(x1 - x2)] + (-1 + 4*t1 + 4*t2 - 8*t1*t2)*
                      Cosh[4*x - 2*(x1 + x2)] + 
                    2*(-1 + 2*t1)*(-1 + 2*t2)*Sinh[4*x - 2*(x1 + x2)])*(-1 - 
                    vb1 - 2*vh1 + (1 + vb1 - 6*vh1)*Cosh[4*(x - xv1)] - 
                    2*(1 + vb1 - 2*vh1)*Sinh[4*(x - xv1)])) + (
              2*(-1 + t1)*Cosh[2*(x - x1)] + (1 - 2*t1)*Sinh[2*(x - x1)])*(
              2*(-1 + t2)*Cosh[2*(x - x2)] + (1 - 2*t2)*
                  Sinh[2*(x - x2)])*(-1 + 3*vb1 + 
                2*vh1 + (1 + 5*vb1 - 2*vh1)*Cosh[4*(x - xv1)] - 
                4*(vb1 - vh1)*Sinh[4*(x - xv1)]) + (
              2*t1*Cosh[2*(x - x1)] + (1 - 2*t1)*Sinh[2*(x - x1)])*(
              2*t2*Cosh[2*(x - x2)] + (1 - 2*t2)*Sinh[2*(x - x2)])*(
              3 - vb1 + 2*vh1 + (5 + vb1 - 2*vh1)*Cosh[4*(x - xv1)] + 
                4*(-1 + vh1)*Sinh[4*(x - xv1)]))]
	
(*  *)
OneOneZOneH[x_, x1_, t1_, x2_, t2_, xv1_, vb1_, vh1_] := Which[
		x<= x1 && x <= x2 && x <= xv1,
		(2*Cosh[2*(x - xv1)]^2*(
            Cosh[2*(x1 - x2)] + (-1 + 4*t1 + 4*t2 - 8*t1*t2)*
                Cosh[4*x - 2*(x1 + x2)] - 
              2*(-1 + 2*t1)*(-1 + 2*t2)*Sinh[4*x - 2*(x1 + x2)])*(
            4*vh1*Tanh[2*(x - xv1)] - (1 + vb1 - 2*vh1)*Tanh[2*(x - xv1)]^2 + 
              2*(2*vh1 + (1 + vb1)*Tanh[2*(-x + xv1)])))/(-((
                  Cosh[2*(x1 - x2)] + (-1 + 4*t1 + 4*t2 - 8*t1*t2)*
                      Cosh[4*x - 2*(x1 + x2)] - 
                    2*(-1 + 2*t1)*(-1 + 2*t2)*Sinh[4*x - 2*(x1 + x2)])*(-1 - 
                    vb1 - 2*vh1 + (1 + vb1 - 6*vh1)*Cosh[4*(x - xv1)] + 
                    2*(1 + vb1 - 2*vh1)*Sinh[4*(x - xv1)])) + (
              2*(-1 + t1)*Cosh[2*(x - x1)] + (-1 + 2*t1)*Sinh[2*(x - x1)])*(
              2*(-1 + t2)*Cosh[2*(x - x2)] + (-1 + 2*t2)*
                  Sinh[2*(x - x2)])*(-1 + 3*vb1 + 
                2*vh1 + (1 + 5*vb1 - 2*vh1)*Cosh[4*(x - xv1)] + 
                4*(vb1 - vh1)*Sinh[4*(x - xv1)]) + (
              2*t1*Cosh[2*(x - x1)] + (-1 + 2*t1)*Sinh[2*(x - x1)])*(
              2*t2*Cosh[2*(x - x2)] + (-1 + 2*t2)*Sinh[2*(x - x2)])*(
              3 - vb1 + 2*vh1 + (5 + vb1 - 2*vh1)*Cosh[4*(x - xv1)] - 
                4*(-1 + vh1)*Sinh[4*(x - xv1)])),
		x<= x1 && x <= x2 && xv1 < x,
		(2*Cosh[2*(x - xv1)]^2*(
            Cosh[2*(x1 - x2)] + (-1 + 4*t1 + 4*t2 - 8*t1*t2)*
                Cosh[4*x - 2*(x1 + x2)] - 
              2*(-1 + 2*t1)*(-1 + 2*t2)*Sinh[4*x - 2*(x1 + x2)])*(
            4*vh1 + 2*(1 + vb1 - 2*vh1)*Tanh[2*(x - xv1)] - (1 + vb1 - 2*vh1)*
                Tanh[2*(x - xv1)]^2))/(-((
                  Cosh[2*(x1 - x2)] + (-1 + 4*t1 + 4*t2 - 8*t1*t2)*
                      Cosh[4*x - 2*(x1 + x2)] - 
                    2*(-1 + 2*t1)*(-1 + 2*t2)*Sinh[4*x - 2*(x1 + x2)])*(-1 - 
                    vb1 - 2*vh1 + (1 + vb1 - 6*vh1)*Cosh[4*(x - xv1)] - 
                    2*(1 + vb1 - 2*vh1)*Sinh[4*(x - xv1)])) + (
              2*(-1 + t1)*Cosh[2*(x - x1)] + (-1 + 2*t1)*Sinh[2*(x - x1)])*(
              2*(-1 + t2)*Cosh[2*(x - x2)] + (-1 + 2*t2)*
                  Sinh[2*(x - x2)])*(-1 + 3*vb1 + 
                2*vh1 + (1 + 5*vb1 - 2*vh1)*Cosh[4*(x - xv1)] - 
                4*(vb1 - vh1)*Sinh[4*(x - xv1)]) + (
              2*t1*Cosh[2*(x - x1)] + (-1 + 2*t1)*Sinh[2*(x - x1)])*(
              2*t2*Cosh[2*(x - x2)] + (-1 + 2*t2)*Sinh[2*(x - x2)])*(
              3 - vb1 + 2*vh1 + (5 + vb1 - 2*vh1)*Cosh[4*(x - xv1)] + 
                4*(-1 + vh1)*Sinh[4*(x - xv1)])),
		x<= x1 && x2 < x && x <= xv1,
		-((((-1 + t1*(4 - 8*t2) + 4*t2)*Cosh[2*(x1 - x2)] + 
                  Cosh[4*x - 2*(x1 + x2)] + 
                  2*(-1 + 2*t1)*(-1 + 2*t2)*Sinh[2*(x1 - x2)])*(-1 - vb1 - 
                  2*vh1 + (1 + vb1 - 6*vh1)*Cosh[4*(x - xv1)] + 
                  2*(1 + vb1 - 2*vh1)*
                    Sinh[4*(x - xv1)]))/(-(((-1 + t1*(4 - 8*t2) + 4*t2)*
                          Cosh[2*(x1 - x2)] + Cosh[4*x - 2*(x1 + x2)] + 
                        2*(-1 + 2*t1)*(-1 + 2*t2)*Sinh[2*(x1 - x2)])*(-1 - 
                        vb1 - 2*vh1 + (1 + vb1 - 6*vh1)*Cosh[4*(x - xv1)] + 
                        2*(1 + vb1 - 2*vh1)*Sinh[4*(x - xv1)])) + (
                  2*(-1 + t1)*Cosh[2*(x - x1)] + (-1 + 2*t1)*
                      Sinh[2*(x - x1)])*(
                  2*(-1 + t2)*Cosh[2*(x - x2)] + (1 - 2*t2)*
                      Sinh[2*(x - x2)])*(-1 + 3*vb1 + 
                    2*vh1 + (1 + 5*vb1 - 2*vh1)*Cosh[4*(x - xv1)] + 
                    4*(vb1 - vh1)*Sinh[4*(x - xv1)]) + (
                  2*t1*Cosh[2*(x - x1)] + (-1 + 2*t1)*Sinh[2*(x - x1)])*(
                  2*t2*Cosh[2*(x - x2)] + (1 - 2*t2)*Sinh[2*(x - x2)])*(
                  3 - vb1 + 2*vh1 + (5 + vb1 - 2*vh1)*Cosh[4*(x - xv1)] - 
                    4*(-1 + vh1)*Sinh[4*(x - xv1)]))),
		x<= x1 && x2 < x && xv1 < x,
		(2*Cosh[2*(x - xv1)]^2*((-1 + t1*(4 - 8*t2) + 4*t2)*Cosh[2*(x1 - x2)] + 
              Cosh[4*x - 2*(x1 + x2)] + 
              2*(-1 + 2*t1)*(-1 + 2*t2)*Sinh[2*(x1 - x2)])*(
            4*vh1 + 2*(1 + vb1 - 2*vh1)*Tanh[2*(x - xv1)] - (1 + vb1 - 2*vh1)*
                Tanh[2*(x - xv1)]^2))/(-(((-1 + t1*(4 - 8*t2) + 4*t2)*
                      Cosh[2*(x1 - x2)] + Cosh[4*x - 2*(x1 + x2)] + 
                    2*(-1 + 2*t1)*(-1 + 2*t2)*Sinh[2*(x1 - x2)])*(-1 - vb1 - 
                    2*vh1 + (1 + vb1 - 6*vh1)*Cosh[4*(x - xv1)] - 
                    2*(1 + vb1 - 2*vh1)*Sinh[4*(x - xv1)])) + (
              2*(-1 + t1)*Cosh[2*(x - x1)] + (-1 + 2*t1)*Sinh[2*(x - x1)])*(
              2*(-1 + t2)*Cosh[2*(x - x2)] + (1 - 2*t2)*
                  Sinh[2*(x - x2)])*(-1 + 3*vb1 + 
                2*vh1 + (1 + 5*vb1 - 2*vh1)*Cosh[4*(x - xv1)] - 
                4*(vb1 - vh1)*Sinh[4*(x - xv1)]) + (
              2*t1*Cosh[2*(x - x1)] + (-1 + 2*t1)*Sinh[2*(x - x1)])*(
              2*t2*Cosh[2*(x - x2)] + (1 - 2*t2)*Sinh[2*(x - x2)])*(
              3 - vb1 + 2*vh1 + (5 + vb1 - 2*vh1)*Cosh[4*(x - xv1)] + 
                4*(-1 + vh1)*Sinh[4*(x - xv1)])),
		x1< x && x <= x2 && x <= xv1,
		(2*Cosh[2*(x - xv1)]^2*((-1 + t1*(4 - 8*t2) + 4*t2)*Cosh[2*(x1 - x2)] + 
              Cosh[4*x - 2*(x1 + x2)] - 
              2*(-1 + 2*t1)*(-1 + 2*t2)*Sinh[2*(x1 - x2)])*(
            4*vh1*Tanh[2*(x - xv1)] - (1 + vb1 - 2*vh1)*Tanh[2*(x - xv1)]^2 + 
              2*(2*vh1 + (1 + vb1)*
                      Tanh[2*(-x + xv1)])))/(-(((-1 + t1*(4 - 8*t2) + 4*t2)*
                      Cosh[2*(x1 - x2)] + Cosh[4*x - 2*(x1 + x2)] - 
                    2*(-1 + 2*t1)*(-1 + 2*t2)*Sinh[2*(x1 - x2)])*(-1 - vb1 - 
                    2*vh1 + (1 + vb1 - 6*vh1)*Cosh[4*(x - xv1)] + 
                    2*(1 + vb1 - 2*vh1)*Sinh[4*(x - xv1)])) + (
              2*(-1 + t1)*Cosh[2*(x - x1)] + (1 - 2*t1)*Sinh[2*(x - x1)])*(
              2*(-1 + t2)*Cosh[2*(x - x2)] + (-1 + 2*t2)*
                  Sinh[2*(x - x2)])*(-1 + 3*vb1 + 
                2*vh1 + (1 + 5*vb1 - 2*vh1)*Cosh[4*(x - xv1)] + 
                4*(vb1 - vh1)*Sinh[4*(x - xv1)]) + (
              2*t1*Cosh[2*(x - x1)] + (1 - 2*t1)*Sinh[2*(x - x1)])*(
              2*t2*Cosh[2*(x - x2)] + (-1 + 2*t2)*Sinh[2*(x - x2)])*(
              3 - vb1 + 2*vh1 + (5 + vb1 - 2*vh1)*Cosh[4*(x - xv1)] - 
                4*(-1 + vh1)*Sinh[4*(x - xv1)])),
		x1< x && x <= x2 && xv1 < x,
		(2*Cosh[2*(x - xv1)]^2*((-1 + t1*(4 - 8*t2) + 4*t2)*Cosh[2*(x1 - x2)] + 
              Cosh[4*x - 2*(x1 + x2)] - 
              2*(-1 + 2*t1)*(-1 + 2*t2)*Sinh[2*(x1 - x2)])*(
            4*vh1 + 2*(1 + vb1 - 2*vh1)*Tanh[2*(x - xv1)] - (1 + vb1 - 2*vh1)*
                Tanh[2*(x - xv1)]^2))/(-(((-1 + t1*(4 - 8*t2) + 4*t2)*
                      Cosh[2*(x1 - x2)] + Cosh[4*x - 2*(x1 + x2)] - 
                    2*(-1 + 2*t1)*(-1 + 2*t2)*Sinh[2*(x1 - x2)])*(-1 - vb1 - 
                    2*vh1 + (1 + vb1 - 6*vh1)*Cosh[4*(x - xv1)] - 
                    2*(1 + vb1 - 2*vh1)*Sinh[4*(x - xv1)])) + (
              2*(-1 + t1)*Cosh[2*(x - x1)] + (1 - 2*t1)*Sinh[2*(x - x1)])*(
              2*(-1 + t2)*Cosh[2*(x - x2)] + (-1 + 2*t2)*
                  Sinh[2*(x - x2)])*(-1 + 3*vb1 + 
                2*vh1 + (1 + 5*vb1 - 2*vh1)*Cosh[4*(x - xv1)] - 
                4*(vb1 - vh1)*Sinh[4*(x - xv1)]) + (
              2*t1*Cosh[2*(x - x1)] + (1 - 2*t1)*Sinh[2*(x - x1)])*(
              2*t2*Cosh[2*(x - x2)] + (-1 + 2*t2)*Sinh[2*(x - x2)])*(
              3 - vb1 + 2*vh1 + (5 + vb1 - 2*vh1)*Cosh[4*(x - xv1)] + 
                4*(-1 + vh1)*Sinh[4*(x - xv1)])),
		x1< x && x2 < x && x <= xv1,
		(2*Cosh[2*(x - xv1)]^2*(
            Cosh[2*(x1 - x2)] + (-1 + 4*t1 + 4*t2 - 8*t1*t2)*
                Cosh[4*x - 2*(x1 + x2)] + 
              2*(-1 + 2*t1)*(-1 + 2*t2)*Sinh[4*x - 2*(x1 + x2)])*(
            4*vh1*Tanh[2*(x - xv1)] - (1 + vb1 - 2*vh1)*Tanh[2*(x - xv1)]^2 + 
              2*(2*vh1 + (1 + vb1)*Tanh[2*(-x + xv1)])))/(-((
                  Cosh[2*(x1 - x2)] + (-1 + 4*t1 + 4*t2 - 8*t1*t2)*
                      Cosh[4*x - 2*(x1 + x2)] + 
                    2*(-1 + 2*t1)*(-1 + 2*t2)*Sinh[4*x - 2*(x1 + x2)])*(-1 - 
                    vb1 - 2*vh1 + (1 + vb1 - 6*vh1)*Cosh[4*(x - xv1)] + 
                    2*(1 + vb1 - 2*vh1)*Sinh[4*(x - xv1)])) + (
              2*(-1 + t1)*Cosh[2*(x - x1)] + (1 - 2*t1)*Sinh[2*(x - x1)])*(
              2*(-1 + t2)*Cosh[2*(x - x2)] + (1 - 2*t2)*
                  Sinh[2*(x - x2)])*(-1 + 3*vb1 + 
                2*vh1 + (1 + 5*vb1 - 2*vh1)*Cosh[4*(x - xv1)] + 
                4*(vb1 - vh1)*Sinh[4*(x - xv1)]) + (
              2*t1*Cosh[2*(x - x1)] + (1 - 2*t1)*Sinh[2*(x - x1)])*(
              2*t2*Cosh[2*(x - x2)] + (1 - 2*t2)*Sinh[2*(x - x2)])*(
              3 - vb1 + 2*vh1 + (5 + vb1 - 2*vh1)*Cosh[4*(x - xv1)] - 
                4*(-1 + vh1)*Sinh[4*(x - xv1)])),
		x1< x && x2 < x && xv1 < x,
		(2*Cosh[2*(x - xv1)]^2*(
            Cosh[2*(x1 - x2)] + (-1 + 4*t1 + 4*t2 - 8*t1*t2)*
                Cosh[4*x - 2*(x1 + x2)] + 
              2*(-1 + 2*t1)*(-1 + 2*t2)*Sinh[4*x - 2*(x1 + x2)])*(
            4*vh1 + 2*(1 + vb1 - 2*vh1)*Tanh[2*(x - xv1)] - (1 + vb1 - 2*vh1)*
                Tanh[2*(x - xv1)]^2))/(-((
                  Cosh[2*(x1 - x2)] + (-1 + 4*t1 + 4*t2 - 8*t1*t2)*
                      Cosh[4*x - 2*(x1 + x2)] + 
                    2*(-1 + 2*t1)*(-1 + 2*t2)*Sinh[4*x - 2*(x1 + x2)])*(-1 - 
                    vb1 - 2*vh1 + (1 + vb1 - 6*vh1)*Cosh[4*(x - xv1)] - 
                    2*(1 + vb1 - 2*vh1)*Sinh[4*(x - xv1)])) + (
              2*(-1 + t1)*Cosh[2*(x - x1)] + (1 - 2*t1)*Sinh[2*(x - x1)])*(
              2*(-1 + t2)*Cosh[2*(x - x2)] + (1 - 2*t2)*
                  Sinh[2*(x - x2)])*(-1 + 3*vb1 + 
                2*vh1 + (1 + 5*vb1 - 2*vh1)*Cosh[4*(x - xv1)] - 
                4*(vb1 - vh1)*Sinh[4*(x - xv1)]) + (
              2*t1*Cosh[2*(x - x1)] + (1 - 2*t1)*Sinh[2*(x - x1)])*(
              2*t2*Cosh[2*(x - x2)] + (1 - 2*t2)*Sinh[2*(x - x2)])*(
              3 - vb1 + 2*vh1 + (5 + vb1 - 2*vh1)*Cosh[4*(x - xv1)] + 
                4*(-1 + vh1)*Sinh[4*(x - xv1)]))]

(*   calculation of chi-list and chi-squares *)
calcchi[pts_List, response_List, guess_List, wp_] :=
	Block[{x1, t1, x2, t2, xv1, vb1, vh1, ptsA, ptsB, ptsH, responseA, 
      responseB, responseH, predA, predB, predH, chilistA, chilistB, 
      chilistH, chisq},
		{responseA, responseB, responseH} = response;
		{ptsA, ptsB, ptsH} = pts;
		{x1, t1, x2, t2, xv1, vb1, vh1} =N[guess,wp];
		predA = N[OneOneZOneA[#1, x1, t1, x2, t2, xv1, vb1, vh1] & /@ ptsA, wp];
		predB = N[OneOneZOneB[#1, x1, t1, x2, t2, xv1, vb1, vh1] & /@ ptsB, wp];
		predH = N[OneOneZOneH[#1, x1, t1, x2, t2, xv1, vb1, vh1] & /@ ptsH, wp];
		chilistA = responseA - predA;
		chilistB = responseB - predB;
		chilistH = responseH - predH;
		chisq =N[ 
        Plus @@ (
            Join[Map[#1^2 &, chilistA], Map[#1^2 &, chilistB], 
              Map[#1^2 &, chilistH]]), wp];
		{{chilistA, chilistB, chilistH}, chisq}
		]

(* ============  MakeOutPutList  ============== *)

MakeOutputList[out_] := 
  Block[
		{x, chisq, param, x1, t1, x2, t2, xv1, vb1, vh1, paramerr, errx1, errt1, 
      errx2, errt2, errxv1, errvb1, errvh1, errorDOF},
		{chisq, param, paramerr, errorDOF, fp} = out;
		{x1, t1, x2, t2, xv1, vb1, vh1} = param;
		If[!(MemberQ[fp, "x1"]),errx1=First[paramerr];
				paramerr=Rest[paramerr], errx1=0];
		If[!(MemberQ[fp, "t1"]),errt1=First[paramerr];paramerr=Rest[paramerr], 
      errt1=0];
		If[!(MemberQ[fp, "x2"]),errx2=First[paramerr];
				paramerr=Rest[paramerr], errx2=0];
		If[!(MemberQ[fp, "t2"]),errt2=First[paramerr];paramerr=Rest[paramerr], 
      errt2=0];
		If[!(MemberQ[fp, "xv1"]),errxv1=First[paramerr];paramerr=Rest[paramerr], 
      errxv1=0];
		If[!(MemberQ[fp, "vb1"]),errvb1=First[paramerr];
				paramerr=Rest[paramerr], errvb1=0];
		If[!(MemberQ[fp, "vh1"]),errvh1=First[paramerr], errvh1=0];
		Print["The model is two gametophyte genes on the male and female gametes \
and a single zygotic viability gene."];
		If[fp=!={},Print["The fixed parameters are ",fp]];
		Print[StringForm[
        "One gametophyte gene position is at `` \[PlusMinus] `` .", 
        NumberForm[x1, 6], NumberForm[errx1, 6]]];
		Print[StringForm[
        "The transmittance of the A gamete by the gene is `` \[PlusMinus] `` \
.", NumberForm[t1, 6], NumberForm[errt1, 6]]];
		Print[StringForm[
        "Another gametophyte gene position is at `` \[PlusMinus] `` .", 
        NumberForm[x2, 6], NumberForm[errx2, 6]]];
		Print[StringForm[
        "The transmittance of the A gamete by the gene is `` \[PlusMinus] `` \
.", NumberForm[t2, 6], NumberForm[errt2, 6]]];
		Print[StringForm[
        "The zygotic viability gene position is at `` \[PlusMinus] `` .", 
        NumberForm[xv1, 6], NumberForm[errxv1, 6]]];
		Print[StringForm[
        "The relative viability of B homozygote is `` \[PlusMinus] `` .", 
        NumberForm[vb1, 6], NumberForm[errvb1, 6]]];
		Print[StringForm[
        "The relative viability of heterozygote is `` \[PlusMinus] `` .", 
        NumberForm[vh1, 6], NumberForm[errvh1, 6]]];
		Print[StringForm[
        "The number of analysed data is `1`.  The root mean square is `2`.", 
        NumberForm[errorDOF + 7-Length[fp], 6], 
        NumberForm[Sqrt[chisq/errorDOF], 6]]]
	]

(*   *)
printprogress[iteration_, chisq_, params_] := 
     Print[
    StringForm[
      "Iteration:`1` ChiSquared:`2`  {x1, t1, x2, t2, xv1, vb1, vh1}:`3`", 
											NumberForm[iteration, 6], NumberForm[chisq, 6], 
      NumberForm[params, 6]]]

(*    *)
DiagonalElements[matrix_] := 
  Module[{k, size = Length[matrix]}, 
   Table[matrix[[k,k]], {k, 1, size}]
		]

(* =========================================== *)
End[]
EndPackage[]