(***********************************************************************
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.
***********************************************************************)

BeginPackage["SegAnalysis`ZTwo`", "SegAnalysis`jacoZ21`", 
  "SegAnalysis`jacoZ22`"]

ZTwoA::usage =
  "ZTwoA[x, xv1, vb1, vh1, xv2, vb2, vh2] gives the expected fraction of A \
genotype at x, when two genes at xv1 and xv2 affect viability of zygote.  The \
relateve viabilities of B homozygote and heterozygote to the A homozygote at \
xv1 and these at xv2 are vb1, vh1, vb2, and vh2, respectively."

ZTwoB::usage =
  "ZTwoB[x, xv1, vb1, vh1, xv2, vb2, vh2] gives the expected fraction of B \
genotype at x, when two genes at xv1 and xv2 affect viability of zygote.  The \
relateve viabilities of B homozygote and heterozygote to the A homozygote at \
xv1 and these at xv2 are vb1, vh1, vb2, and vh2, respectively."

ZTwoH::usage =
  "ZTwoH[x, xv1, vb1, vh1, xv2, vb2, vh2] gives the expected fraction of H \
genotype at x, when two genes at xv1 and xv2 affect viability of zygote.  The \
relateve viabilities of B homozygote and heterozygote to the A homozygote at \
xv1 and these at xv2 are vb1, vh1, vb2, and vh2, respectively."

Begin["`Private`"]

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

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

(* =========    ZTwo     =========== *)
(* This uses Guess-Newton method for finding the least-
    squares solution for two zygotic viability genes model  *)

mnlr[nptsA_, nptsB_, nptsH_, nresponseA_, nresponseB_, nresponseH_, xv1_, 
    vb1_, vh1_, xv2_, vb2_, vh2_, fp_, maxits_, wp_, pg_, ag_, mdls_, 
    zt_] := 
  Module[{predA, predB, predH, chilistA, chilistB, chilistH, tmpchilistA, 
      tmpchilistB, tmpchilistH, tmpchilist, tmpnptsA, tmpnptsB, tmpnptsH, 
      chisq, derxv1,dervb1, dervh1, derxv2, dervb2, dervh2,der, alpha, beta, 
      chisq0, newchi, tmpchi, its, accDelta, NAcc, precDelta, NPrec, 
      invaliddigits, accflag0, accflag1 = False, accflag2 = False, precflag0, 
      precflag1 = False, precflag2 = False, ovflag = False, guessxv1 = xv1,
      guessvb1 = vb1, guessvh1 = vh1, guessxv2 = xv2,guessvb2 = vb2, 
      guessvh2 = vh2, oldxv1, oldvb1, oldvh1, oldxv2, oldvb2, oldvh2, delta, 
      deltaxv1, deltavb1, deltavh1, deltaxv2, deltavb2, deltavh2, tmpxv1, 
      tmpvb1, tmpvh1, tmpxv2, tmpvb2, tmpvh2, 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}, {
          guessxv1,guessvb1, guessvh1, guessxv2,guessvb2, guessvh2}, wp];
		tmpchi = chisq; newchi = -1; its = 1; 
		oldxv1 = guessxv1;
		oldvb1 = guessvb1; oldvh1 = guessvh1; 
		oldxv2 = guessxv2;
		oldvb2 = guessvb2; oldvh2 = guessvh2;

		(* 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;
					oldxv1 = guessxv1; 
					oldvb1 = guessvb1; oldvh1 = guessvh1;
					oldxv2 = guessxv2; 
					oldvb2 = guessvb2; oldvh2 = guessvh2;
					
					(* === find deltaxv1, deltava1, 
            deltavb1 and deltavh1 by LinearSolve normalmatrix ==== *)
					
					(* === delete not differentiable position ==== *)
					tmpchilistA = Delete[chilistA, Position[nptsA, guessxv1]];
					tmpnptsA =Delete[nptsA, Position[nptsA, guessxv1]];
					tmpchilistB = Delete[chilistB, Position[nptsB, guessxv1]];
					tmpnptsB =Delete[nptsB, Position[nptsB, guessxv1]];
					tmpchilistH = Delete[chilistH, Position[nptsH, guessxv1]];
					tmpnptsH =Delete[nptsH, Position[nptsH, guessxv1]];
					tmpchilistA = Delete[tmpchilistA, Position[tmpnptsA, guessxv2]];
					tmpnptsA =Delete[tmpnptsA, Position[tmpnptsA, guessxv2]];
					tmpchilistB = Delete[tmpchilistB, Position[tmpnptsB, guessxv2]];
					tmpnptsB =Delete[tmpnptsB, Position[tmpnptsB, guessxv2]];
					tmpchilistH = Delete[tmpchilistH, Position[tmpnptsH, guessxv2]];
					tmpnptsH =Delete[tmpnptsH, Position[tmpnptsH, guessxv2]];
					tmpchilist = Join[tmpchilistA , tmpchilistB, tmpchilistH];
					(* calculate jacobian *)
					der={};
					delta={};
					
					If[!(MemberQ[fp, "xv1"]),
					derxv1 =
              SegAnalysis`jacoZ21`Private`a[guessxv1, guessvb1, guessvh1, 
                guessxv2, guessvb2, guessvh2, tmpnptsA, tmpnptsB, tmpnptsH, 
                wp];
							der=Append[der,derxv1]];
					
					If[!(MemberQ[fp, "vb1"]),
						dervb1 =
              SegAnalysis`jacoZ21`Private`b[guessxv1, guessvb1, guessvh1, 
                guessxv2, guessvb2, guessvh2, tmpnptsA, tmpnptsB, tmpnptsH, 
                wp];
						der=Append[der,dervb1]];
					
					If[!(MemberQ[fp, "vh1"]),
            dervh1 =SegAnalysis`jacoZ21`Private`c[guessxv1, guessvb1, 
                guessvh1, guessxv2, guessvb2, guessvh2, tmpnptsA, tmpnptsB, 
                tmpnptsH, wp];
						der=Append[der,dervh1]];
					
					If[!(MemberQ[fp, "xv2"]),
					derxv2 =
              SegAnalysis`jacoZ22`Private`a[guessxv1, guessvb1, guessvh1, 
                guessxv2, guessvb2, guessvh2, tmpnptsA, tmpnptsB, tmpnptsH, 
                wp];
							der=Append[der,derxv2]];
					
					If[!(MemberQ[fp, "vb2"]),
						dervb2 =
              SegAnalysis`jacoZ22`Private`b[guessxv1, guessvb1, guessvh1, 
                guessxv2, guessvb2, guessvh2, tmpnptsA, tmpnptsB, tmpnptsH, 
                wp];
						der=Append[der,dervb2]];
					
					If[!(MemberQ[fp, "vh2"]),
            dervh2 =SegAnalysis`jacoZ22`Private`c[guessxv1, guessvb1, 
                guessvh1, guessxv2, guessvb2, guessvh2, tmpnptsA, tmpnptsB, 
                tmpnptsH, wp];
						der=Append[der,dervh2]];
									
					
					(* 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, "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];delta=Rest[delta],  
            deltavh1=0];
					If[!(MemberQ[fp, "xv2"]),deltaxv2=First[delta];delta=Rest[delta], 
            deltaxv2=0];
					If[!(MemberQ[fp, "vb2"]),deltavb2=First[delta]; delta=Rest[delta],  
            deltavb2=0];
					If[!(MemberQ[fp, "vh2"]),deltavh2=First[delta],  deltavh2=0];
					
					{tmpxv1,tmpvb1, tmpvh1, tmpxv2,tmpvb2, tmpvh2} = 
            N[{oldxv1, oldvb1, oldvh1, oldxv2, oldvb2, oldvh2} +  {deltaxv1, 
                  deltavb1, deltavh1, deltaxv2, deltavb2, deltavh2}, wp];
					
					{{newchilistA, newchilistB, newchilistH}, newchi} = 
            calcchi[{nptsA, nptsB, nptsH}, {nresponseA, nresponseB, 
                nresponseH}, {tmpxv1,tmpvb1, tmpvh1, tmpxv2,tmpvb2, tmpvh2}, 
              wp];
					
					(* find newguess for smaller chisq  *)
					dits = 1; dmaxits = 15;
					While[newchi - chisq >= 0 && dits <= dmaxits,
						{tmpxv1,tmpvb1, tmpvh1, tmpxv2,tmpvb2, tmpvh2} = 
              N[{oldxv1, oldvb1, oldvh1, oldxv2, oldvb2, 
                    oldvh2} + (2^(-1*dits))*{deltaxv1, deltavb1, deltavh1, 
                      deltaxv2, deltavb2, deltavh2}, wp];
						{{newchilistA, newchilistB, newchilistH}, newchi} = 
              calcchi[{nptsA, nptsB, nptsH}, {nresponseA, nresponseB, 
                  nresponseH}, {tmpxv1,tmpvb1, tmpvh1, tmpxv2,tmpvb2, 
                  tmpvh2}, wp];
						dits++];
					If[dits > dmaxits, Message[SegregAnalysis::lmincrease, dmaxits]];
					tmpchi = newchi;
					{guessxv1, guessvb1, guessvh1, guessxv2, guessvb2, guessvh2} = {tmpxv1,
              tmpvb1, tmpvh1, tmpxv2,tmpvb2, tmpvh2};
					chilistA = newchilistA;
					chilistB = newchilistB;
					chilistH = newchilistH;
					printprogress[its, 
            newchi, {guessxv1, guessvb1, guessvh1, guessxv2, guessvb2, 
              guessvh2}];
					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] - 6+Length[fp];
		paramerr = Sqrt[ ( newchi/errorDOF) DiagonalElements[Inverse[alpha]]];
		{newchi, {guessxv1, guessvb1, guessvh1, guessxv2, guessvb2, guessvh2}, 
      paramerr, errorDOF, fp}
	]  (* end ZTwo *)

(*   *)
ZTwoA[x_, xv1_, vb1_, vh1_, xv2_, vb2_, vh2_] := Which[
		x <= xv1,
		(4*Tanh[2*(x - xv1)]*(
              4 + 2*(2 + vh1 + vb2*vh1 - 2*vh2 - 2*vh1*vh2)*
                  Tanh[2*(xv1 - xv2)] + (1 + vb2)*(1 + vh1)*
                  Tanh[2*(xv1 - xv2)]^2) + 
          Tanh[2*(x - xv1)]^2*(
              4*(1 + vb1*vb2 - 2*vh1*vh2) + 
                4*(1 + vb1*vb2 + vh1 + vb2*vh1)*
                  Tanh[2*(xv1 - xv2)] + (1 + vb1 + 2*vh1)*(1 + vb2 - 2*vh2)*
                  Tanh[2*(xv1 - xv2)]^2) + 
          4*(4 + 4*Tanh[2*xv1 - 2*xv2] - (1 + vb1 + 2*vh1)*vh2*
                  Tanh[2*(x - xv1)]^2*Tanh[2*(xv1 - xv2)] + 
                Tanh[2*(xv1 - xv2)]^2 + vb2*Tanh[2*(xv1 - xv2)]^2 - 
                2*vh2*Tanh[2*(xv1 - xv2)]^2 + 4*vh2*Tanh[2*(-xv1 + xv2)] + 
                2*vh2*Tanh[2*(-x + xv1)]*(
                    2*vh1 + vh1*Tanh[2*(xv1 - xv2)]^2 + 
                      Tanh[2*(-xv1 + xv2)]^2)))/(
        4*(4*(1 + vb1*vb2 + 2*vh1*vh2) + 
              4*(1 + vb1*vb2 + 2*vh1*vh2)*
                Tanh[2*(xv1 - xv2)] + (
                  vb2 - 2*vh1 - 2*vb2*vh1 + vb1*(1 + vb2 - 2*vh2) - 2*vh2 + 
                    4*vh1*vh2)*Tanh[2*(xv1 - xv2)]^2 + 
              4*(vh1 + vb2*vh1 + vh2 + vb1*vh2)*Tanh[2*(-xv1 + xv2)] + 
              Tanh[2*(-xv1 + xv2)]^2)),
		xv1 < x && x <= xv2,
		((3 - vb1 + 2*vh1 + (5 + vb1 - 2*vh1)*Cosh[4*(x - xv1)] + 
              4*(-1 + vh1)*Sinh[4*(x - xv1)])*(
            3 - vb2 + 2*vh2 + (5 + vb2 - 2*vh2)*Cosh[4*(x - xv2)] - 
              4*(-1 + vh2)*Sinh[4*(x - xv2)]))/(
        2*(6 - 2*vb1 - 2*vb2 + 6*vb1*vb2 + 4*vh1 + 4*vb2*vh1 + 4*vh2 + 
              4*vb1*vh2 + 8*vh1*vh2 + 
              2*(3 - vb1 - vb2 + 3*vb1*vb2 + 2*vh1 + 2*vb2*vh1 + 2*vh2 + 
                    2*vb1*vh2 + 4*vh1*vh2)*Cosh[4*(x - xv1)] + 
              2*(3 - vb1 - vb2 + 3*vb1*vb2 + 2*vh1 + 2*vb2*vh1 + 2*vh2 + 
                    2*vb1*vh2 + 4*vh1*vh2)*Cosh[4*(x - xv2)] + 
              13*Cosh[4*(xv1 - xv2)] + 5*vb1*Cosh[4*(xv1 - xv2)] + 
              5*vb2*Cosh[4*(xv1 - xv2)] + 13*vb1*vb2*Cosh[4*(xv1 - xv2)] - 
              14*vh1*Cosh[4*(xv1 - xv2)] - 14*vb2*vh1*Cosh[4*(xv1 - xv2)] - 
              14*vh2*Cosh[4*(xv1 - xv2)] - 14*vb1*vh2*Cosh[4*(xv1 - xv2)] + 
              36*vh1*vh2*Cosh[4*(xv1 - xv2)] + Cosh[8*x - 4*(xv1 + xv2)] + 
              vb1*Cosh[8*x - 4*(xv1 + xv2)] + vb2*Cosh[8*x - 4*(xv1 + xv2)] + 
              vb1*vb2*Cosh[8*x - 4*(xv1 + xv2)] + 
              2*vh1*Cosh[8*x - 4*(xv1 + xv2)] + 
              2*vb2*vh1*Cosh[8*x - 4*(xv1 + xv2)] + 
              2*vh2*Cosh[8*x - 4*(xv1 + xv2)] + 
              2*vb1*vh2*Cosh[8*x - 4*(xv1 + xv2)] + 
              4*vh1*vh2*Cosh[8*x - 4*(xv1 + xv2)] - 4*Sinh[4*(x - xv1)] + 
              4*vb1*Sinh[4*(x - xv1)] + 4*vb2*Sinh[4*(x - xv1)] - 
              4*vb1*vb2*Sinh[4*(x - xv1)] + 4*Sinh[4*(x - xv2)] - 
              4*vb1*Sinh[4*(x - xv2)] - 4*vb2*Sinh[4*(x - xv2)] + 
              4*vb1*vb2*Sinh[4*(x - xv2)] + 12*Sinh[4*(xv1 - xv2)] + 
              4*vb1*Sinh[4*(xv1 - xv2)] + 4*vb2*Sinh[4*(xv1 - xv2)] + 
              12*vb1*vb2*Sinh[4*(xv1 - xv2)] - 16*vh1*Sinh[4*(xv1 - xv2)] - 
              16*vb2*vh1*Sinh[4*(xv1 - xv2)] - 16*vh2*Sinh[4*(xv1 - xv2)] - 
              16*vb1*vh2*Sinh[4*(xv1 - xv2)] + 
              32*vh1*vh2*Sinh[4*(xv1 - xv2)])),
		xv2 < x,
		(4*(4 + (1 + vb1*vb2 + vh2 + vb1*vh2)*Tanh[2*(x - xv2)]^2)*
            Tanh[2*(xv1 - xv2)] + 
          4*(4 + Tanh[2*(xv1 - xv2)]^2 + vb1*Tanh[2*(xv1 - xv2)]^2 - 
                2*vh1*Tanh[2*(xv1 - xv2)]^2 + 4*vh1*Tanh[2*(-xv1 + xv2)]) + 
          Tanh[2*(x - xv2)]^2*(
              4 + 4*vb1*vb2 - 8*vh1*vh2 + Tanh[2*(xv1 - xv2)]^2 + 
                vb1*Tanh[2*(xv1 - xv2)]^2 + vb2*Tanh[2*(xv1 - xv2)]^2 + 
                vb1*vb2*Tanh[2*(xv1 - xv2)]^2 - 2*vh1*Tanh[2*(xv1 - xv2)]^2 - 
                2*vb2*vh1*Tanh[2*(xv1 - xv2)]^2 + 
                2*vh2*Tanh[2*(xv1 - xv2)]^2 + 
                2*vb1*vh2*Tanh[2*(xv1 - xv2)]^2 - 
                4*vh1*vh2*Tanh[2*(xv1 - xv2)]^2 + 
                4*vh1*Tanh[2*(-xv1 + xv2)] + 4*vb2*vh1*Tanh[2*(-xv1 + xv2)] + 
                8*vh1*vh2*Tanh[2*(-xv1 + xv2)]) + 
          4*Tanh[2*(x - xv2)]*(
              4*vh1*(1 + vh2)*
                  Tanh[2*(xv1 - xv2)] - (1 + vb1 - 2*vh1)*(1 + vh2)*
                  Tanh[2*(xv1 - xv2)]^2 + 
                2*(-2 + 2*vh1*vh2 + 2*Tanh[2*(-xv1 + xv2)] + 
                      vh2*Tanh[2*(-xv1 + xv2)] + 
                      vb1*vh2*Tanh[2*(-xv1 + xv2)])))/(
        4*(4*(1 + vb1*vb2 + 2*vh1*vh2) + 
              4*(1 + vb1*vb2 + 2*vh1*vh2)*
                Tanh[2*(xv1 - xv2)] + (
                  vb2 - 2*vh1 - 2*vb2*vh1 + vb1*(1 + vb2 - 2*vh2) - 2*vh2 + 
                    4*vh1*vh2)*Tanh[2*(xv1 - xv2)]^2 + 
              4*(vh1 + vb2*vh1 + vh2 + vb1*vh2)*Tanh[2*(-xv1 + xv2)] + 
              Tanh[2*(-xv1 + xv2)]^2))]
	
(*  *)
ZTwoB[x_, xv1_, vb1_, vh1_, xv2_, vb2_, vh2_] := Which[
		x <= xv1,
		(4*Tanh[2*(x - xv1)]*(
              4*vb1*vb2 + 
                2*(2*vb1*vb2 + vh1 + vb2*vh1 - 2*vb1*vh2 - 2*vh1*vh2)*
                  Tanh[2*(xv1 - xv2)] + (1 + vb2)*(vb1 + vh1)*
                  Tanh[2*(xv1 - xv2)]^2) + 
          Tanh[2*(x - xv1)]^2*(
              4*(1 + vb1*vb2 - 2*vh1*vh2) + 
                4*(1 + vb1*vb2 + vh1 + vb2*vh1)*
                  Tanh[2*(xv1 - xv2)] + (1 + vb1 + 2*vh1)*(1 + vb2 - 2*vh2)*
                  Tanh[2*(xv1 - xv2)]^2) + 
          4*(-((1 + vb1 + 2*vh1)*vh2*Tanh[2*(x - xv1)]^2*
                      Tanh[2*(xv1 - xv2)]) + 
                2*vh2*Tanh[2*(-x + xv1)]*(
                    2*vh1 + vb1*Tanh[2*(xv1 - xv2)]^2 + 
                      vh1*Tanh[2*(xv1 - xv2)]^2) + 
                vb1*(4*vb2 + 4*vb2*Tanh[2*xv1 - 2*xv2] + 
                      Tanh[2*(xv1 - xv2)]^2 + vb2*Tanh[2*(xv1 - xv2)]^2 - 
                      2*vh2*Tanh[2*(xv1 - xv2)]^2 + 
                      4*vh2*Tanh[2*(-xv1 + xv2)])))/(
        4*(4*(1 + vb1*vb2 + 2*vh1*vh2) + 
              4*(1 + vb1*vb2 + 2*vh1*vh2)*
                Tanh[2*(xv1 - xv2)] + (
                  vb2 - 2*vh1 - 2*vb2*vh1 + vb1*(1 + vb2 - 2*vh2) - 2*vh2 + 
                    4*vh1*vh2)*Tanh[2*(xv1 - xv2)]^2 + 
              4*(vh1 + vb2*vh1 + vh2 + vb1*vh2)*Tanh[2*(-xv1 + xv2)] + 
              Tanh[2*(-xv1 + xv2)]^2)),
		xv1 < x && x <= xv2,
		((-1 + 3*vb1 + 2*vh1 + (1 + 5*vb1 - 2*vh1)*Cosh[4*(x - xv1)] - 
              4*(vb1 - vh1)*Sinh[4*(x - xv1)])*(-1 + 3*vb2 + 
              2*vh2 + (1 + 5*vb2 - 2*vh2)*Cosh[4*(x - xv2)] + 
              4*(vb2 - vh2)*Sinh[4*(x - xv2)]))/(
        2*(6 - 2*vb1 - 2*vb2 + 6*vb1*vb2 + 4*vh1 + 4*vb2*vh1 + 4*vh2 + 
              4*vb1*vh2 + 8*vh1*vh2 + 
              2*(3 - vb1 - vb2 + 3*vb1*vb2 + 2*vh1 + 2*vb2*vh1 + 2*vh2 + 
                    2*vb1*vh2 + 4*vh1*vh2)*Cosh[4*(x - xv1)] + 
              2*(3 - vb1 - vb2 + 3*vb1*vb2 + 2*vh1 + 2*vb2*vh1 + 2*vh2 + 
                    2*vb1*vh2 + 4*vh1*vh2)*Cosh[4*(x - xv2)] + 
              13*Cosh[4*(xv1 - xv2)] + 5*vb1*Cosh[4*(xv1 - xv2)] + 
              5*vb2*Cosh[4*(xv1 - xv2)] + 13*vb1*vb2*Cosh[4*(xv1 - xv2)] - 
              14*vh1*Cosh[4*(xv1 - xv2)] - 14*vb2*vh1*Cosh[4*(xv1 - xv2)] - 
              14*vh2*Cosh[4*(xv1 - xv2)] - 14*vb1*vh2*Cosh[4*(xv1 - xv2)] + 
              36*vh1*vh2*Cosh[4*(xv1 - xv2)] + Cosh[8*x - 4*(xv1 + xv2)] + 
              vb1*Cosh[8*x - 4*(xv1 + xv2)] + vb2*Cosh[8*x - 4*(xv1 + xv2)] + 
              vb1*vb2*Cosh[8*x - 4*(xv1 + xv2)] + 
              2*vh1*Cosh[8*x - 4*(xv1 + xv2)] + 
              2*vb2*vh1*Cosh[8*x - 4*(xv1 + xv2)] + 
              2*vh2*Cosh[8*x - 4*(xv1 + xv2)] + 
              2*vb1*vh2*Cosh[8*x - 4*(xv1 + xv2)] + 
              4*vh1*vh2*Cosh[8*x - 4*(xv1 + xv2)] - 4*Sinh[4*(x - xv1)] + 
              4*vb1*Sinh[4*(x - xv1)] + 4*vb2*Sinh[4*(x - xv1)] - 
              4*vb1*vb2*Sinh[4*(x - xv1)] + 4*Sinh[4*(x - xv2)] - 
              4*vb1*Sinh[4*(x - xv2)] - 4*vb2*Sinh[4*(x - xv2)] + 
              4*vb1*vb2*Sinh[4*(x - xv2)] + 12*Sinh[4*(xv1 - xv2)] + 
              4*vb1*Sinh[4*(xv1 - xv2)] + 4*vb2*Sinh[4*(xv1 - xv2)] + 
              12*vb1*vb2*Sinh[4*(xv1 - xv2)] - 16*vh1*Sinh[4*(xv1 - xv2)] - 
              16*vb2*vh1*Sinh[4*(xv1 - xv2)] - 16*vh2*Sinh[4*(xv1 - xv2)] - 
              16*vb1*vh2*Sinh[4*(xv1 - xv2)] + 
              32*vh1*vh2*Sinh[4*(xv1 - xv2)])),
		xv2 < x,
		(-4*Tanh[2*(x - xv2)]*(
              4*vb1*vb2 - 4*vh1*vh2 + 
                2*(2*vb1*vb2 - 2*vb2*vh1 + vh2 + vb1*vh2 - 2*vh1*vh2)*
                  Tanh[2*(xv1 - xv2)] + (1 + vb1 - 2*vh1)*(vb2 + vh2)*
                  Tanh[2*(xv1 - xv2)]^2) + 
          4*vb2*(4*vb1 + 4*vb1*Tanh[2*xv1 - 2*xv2] + Tanh[2*(xv1 - xv2)]^2 + 
                vb1*Tanh[2*(xv1 - xv2)]^2 - 2*vh1*Tanh[2*(xv1 - xv2)]^2 + 
                4*vh1*Tanh[2*(-xv1 + xv2)]) + 
          Tanh[2*(x - xv2)]^2*(
              4 + 4*vb1*vb2 - 8*vh1*vh2 + 4*Tanh[2*xv1 - 2*xv2] + 
                4*vb1*vb2*Tanh[2*xv1 - 2*xv2] + 4*vh2*Tanh[2*xv1 - 2*xv2] + 
                4*vb1*vh2*Tanh[2*xv1 - 2*xv2] + Tanh[2*(xv1 - xv2)]^2 + 
                vb1*Tanh[2*(xv1 - xv2)]^2 + vb2*Tanh[2*(xv1 - xv2)]^2 + 
                vb1*vb2*Tanh[2*(xv1 - xv2)]^2 - 2*vh1*Tanh[2*(xv1 - xv2)]^2 - 
                2*vb2*vh1*Tanh[2*(xv1 - xv2)]^2 + 
                2*vh2*Tanh[2*(xv1 - xv2)]^2 + 
                2*vb1*vh2*Tanh[2*(xv1 - xv2)]^2 - 
                4*vh1*vh2*Tanh[2*(xv1 - xv2)]^2 + 
                4*vh1*Tanh[2*(-xv1 + xv2)] + 4*vb2*vh1*Tanh[2*(-xv1 + xv2)] + 
                8*vh1*vh2*Tanh[2*(-xv1 + xv2)]))/(
        4*(4*(1 + vb1*vb2 + 2*vh1*vh2) + 
              4*(1 + vb1*vb2 + 2*vh1*vh2)*
                Tanh[2*(xv1 - xv2)] + (
                  vb2 - 2*vh1 - 2*vb2*vh1 + vb1*(1 + vb2 - 2*vh2) - 2*vh2 + 
                    4*vh1*vh2)*Tanh[2*(xv1 - xv2)]^2 + 
              4*(vh1 + vb2*vh1 + vh2 + vb1*vh2)*Tanh[2*(-xv1 + xv2)] + 
              Tanh[2*(-xv1 + xv2)]^2))]
	
(*  *)
ZTwoH[x_, xv1_, vb1_, vh1_, xv2_, vb2_, vh2_] := Which[
		x <= xv1,
		(-(Tanh[2*(x - xv1)]^2*(
                  4*(1 + vb1*vb2 - 2*vh1*vh2) - 
                    4*(1 + vb1 + 2*vh1)*vh2*
                      Tanh[2*(xv1 - xv2)] + (1 + vb1 + 2*vh1)*(
                        1 + vb2 - 2*vh2)*Tanh[2*(xv1 - xv2)]^2)) + 
          4*Tanh[2*(x - xv1)]*(
              4*vh1*vh2 - 
                2*(1 + vb1*vb2 + vh1 + vb2*vh1 - vh2 - vb1*vh2 - 2*vh1*vh2)*
                  Tanh[2*(xv1 - xv2)] + (1 + vb1 + 2*vh1)*vh2*
                  Tanh[2*(xv1 - xv2)]^2) + 
          2*(-2*(1 + vb1*vb2 + vh1 + vb2*vh1)*Tanh[2*(x - xv1)]^2*
                  Tanh[2*(xv1 - xv2)] + 
                2*vh1*(4*vh2 + 4*vh2*Tanh[2*xv1 - 2*xv2] - 
                      Tanh[2*(xv1 - xv2)]^2 - vb2*Tanh[2*(xv1 - xv2)]^2 + 
                      2*vh2*Tanh[2*(xv1 - xv2)]^2 + 2*Tanh[2*(-xv1 + xv2)] + 
                      2*vb2*Tanh[2*(-xv1 + xv2)]) + 
                Tanh[2*(-x + xv1)]*(
                    4 + 4*vb1*vb2 + vb1*Tanh[2*(xv1 - xv2)]^2 + 
                      vb2*Tanh[2*(xv1 - xv2)]^2 + 
                      vb1*vb2*Tanh[2*(xv1 - xv2)]^2 + 
                      2*vh1*Tanh[2*(xv1 - xv2)]^2 + 
                      2*vb2*vh1*Tanh[2*(xv1 - xv2)]^2 + 
                      Tanh[2*(-xv1 + xv2)]^2)))/(
        2*(4*(1 + vb1*vb2 + 2*vh1*vh2) + 
              4*(1 + vb1*vb2 + 2*vh1*vh2)*
                Tanh[2*(xv1 - xv2)] + (
                  vb2 - 2*vh1 - 2*vb2*vh1 + vb1*(1 + vb2 - 2*vh2) - 2*vh2 + 
                    4*vh1*vh2)*Tanh[2*(xv1 - xv2)]^2 + 
              4*(vh1 + vb2*vh1 + vh2 + vb1*vh2)*Tanh[2*(-xv1 + xv2)] + 
              Tanh[2*(-xv1 + xv2)]^2)),
		xv1 < x && x <= xv2,
		((1 + vb1 + 2*vh1 - (1 + vb1 - 6*vh1)*Cosh[4*(x - xv1)] + 
              2*(1 + vb1 - 2*vh1)*Sinh[4*(x - xv1)])*(
            1 + vb2 + 2*vh2 - (1 + vb2 - 6*vh2)*Cosh[4*(x - xv2)] - 
              2*(1 + vb2 - 2*vh2)*Sinh[4*(x - xv2)]))/(
        6 - 2*vb1 - 2*vb2 + 6*vb1*vb2 + 4*vh1 + 4*vb2*vh1 + 4*vh2 + 
          4*vb1*vh2 + 8*vh1*vh2 + 
          2*(3 - vb1 - vb2 + 3*vb1*vb2 + 2*vh1 + 2*vb2*vh1 + 2*vh2 + 
                2*vb1*vh2 + 4*vh1*vh2)*Cosh[4*(x - xv1)] + 
          2*(3 - vb1 - vb2 + 3*vb1*vb2 + 2*vh1 + 2*vb2*vh1 + 2*vh2 + 
                2*vb1*vh2 + 4*vh1*vh2)*Cosh[4*(x - xv2)] + 
          13*Cosh[4*(xv1 - xv2)] + 5*vb1*Cosh[4*(xv1 - xv2)] + 
          5*vb2*Cosh[4*(xv1 - xv2)] + 13*vb1*vb2*Cosh[4*(xv1 - xv2)] - 
          14*vh1*Cosh[4*(xv1 - xv2)] - 14*vb2*vh1*Cosh[4*(xv1 - xv2)] - 
          14*vh2*Cosh[4*(xv1 - xv2)] - 14*vb1*vh2*Cosh[4*(xv1 - xv2)] + 
          36*vh1*vh2*Cosh[4*(xv1 - xv2)] + Cosh[8*x - 4*(xv1 + xv2)] + 
          vb1*Cosh[8*x - 4*(xv1 + xv2)] + vb2*Cosh[8*x - 4*(xv1 + xv2)] + 
          vb1*vb2*Cosh[8*x - 4*(xv1 + xv2)] + 
          2*vh1*Cosh[8*x - 4*(xv1 + xv2)] + 
          2*vb2*vh1*Cosh[8*x - 4*(xv1 + xv2)] + 
          2*vh2*Cosh[8*x - 4*(xv1 + xv2)] + 
          2*vb1*vh2*Cosh[8*x - 4*(xv1 + xv2)] + 
          4*vh1*vh2*Cosh[8*x - 4*(xv1 + xv2)] - 4*Sinh[4*(x - xv1)] + 
          4*vb1*Sinh[4*(x - xv1)] + 4*vb2*Sinh[4*(x - xv1)] - 
          4*vb1*vb2*Sinh[4*(x - xv1)] + 4*Sinh[4*(x - xv2)] - 
          4*vb1*Sinh[4*(x - xv2)] - 4*vb2*Sinh[4*(x - xv2)] + 
          4*vb1*vb2*Sinh[4*(x - xv2)] + 12*Sinh[4*(xv1 - xv2)] + 
          4*vb1*Sinh[4*(xv1 - xv2)] + 4*vb2*Sinh[4*(xv1 - xv2)] + 
          12*vb1*vb2*Sinh[4*(xv1 - xv2)] - 16*vh1*Sinh[4*(xv1 - xv2)] - 
          16*vb2*vh1*Sinh[4*(xv1 - xv2)] - 16*vh2*Sinh[4*(xv1 - xv2)] - 
          16*vb1*vh2*Sinh[4*(xv1 - xv2)] + 32*vh1*vh2*Sinh[4*(xv1 - xv2)]),
		xv2 < x,
		(-(Tanh[2*(x - xv2)]^2*(
                  4*(1 + vb1*vb2 - 2*vh1*vh2) + 
                    4*(1 + vb1*vb2 - vh1 - vb2*vh1 + vh2 + vb1*vh2 - 
                          2*vh1*vh2)*
                      Tanh[2*(xv1 - xv2)] + (1 + vb1 - 2*vh1)*(
                        1 + vb2 + 2*vh2)*Tanh[2*(xv1 - xv2)]^2)) + 
          4*vh2*(4*vh1 + 4*vh1*Tanh[2*xv1 - 2*xv2] - Tanh[2*(xv1 - xv2)]^2 - 
                vb1*Tanh[2*(xv1 - xv2)]^2 + 2*vh1*Tanh[2*(xv1 - xv2)]^2 + 
                2*Tanh[2*(-xv1 + xv2)] + 2*vb1*Tanh[2*(-xv1 + xv2)]) + 
          2*Tanh[2*(x - xv2)]*(
              4 + 4*vb1*vb2 - 8*vh1*vh2 + 4*Tanh[2*(xv1 - xv2)] + 
                4*vb1*vb2*Tanh[2*(xv1 - xv2)] + 4*vh2*Tanh[2*(xv1 - xv2)] + 
                4*vb1*vh2*Tanh[2*(xv1 - xv2)] + Tanh[2*(xv1 - xv2)]^2 + 
                vb1*Tanh[2*(xv1 - xv2)]^2 + vb2*Tanh[2*(xv1 - xv2)]^2 + 
                vb1*vb2*Tanh[2*(xv1 - xv2)]^2 - 2*vh1*Tanh[2*(xv1 - xv2)]^2 - 
                2*vb2*vh1*Tanh[2*(xv1 - xv2)]^2 + 
                2*vh2*Tanh[2*(xv1 - xv2)]^2 + 
                2*vb1*vh2*Tanh[2*(xv1 - xv2)]^2 - 
                4*vh1*vh2*Tanh[2*(xv1 - xv2)]^2 + 
                4*vh1*Tanh[2*(-xv1 + xv2)] + 4*vb2*vh1*Tanh[2*(-xv1 + xv2)] + 
                8*vh1*vh2*Tanh[2*(-xv1 + xv2)]))/(
        2*(4*(1 + vb1*vb2 + 2*vh1*vh2) + 
              4*(1 + vb1*vb2 + 2*vh1*vh2)*
                Tanh[2*(xv1 - xv2)] + (
                  vb2 - 2*vh1 - 2*vb2*vh1 + vb1*(1 + vb2 - 2*vh2) - 2*vh2 + 
                    4*vh1*vh2)*Tanh[2*(xv1 - xv2)]^2 + 
              4*(vh1 + vb2*vh1 + vh2 + vb1*vh2)*Tanh[2*(-xv1 + xv2)] + 
              Tanh[2*(-xv1 + xv2)]^2))]

(*   calculation of chi-list and chi-squares *)
calcchi[pts_List, response_List, guess_List, wp_] :=
	Block[{xv1, vb1, vh1, xv2, vb2, vh2, ptsA, ptsB, ptsH, responseA, responseB, 
      responseH, predA, predB, predH, chilistA, chilistB, chilistH, chisq},
		{responseA, responseB, responseH} = response;
		{ptsA, ptsB, ptsH} = pts;
		{xv1, vb1, vh1, xv2, vb2, vh2} =N[guess,wp];
		predA = N[ZTwoA[#1, xv1, vb1, vh1, xv2, vb2, vh2] & /@ ptsA, wp];
		predB = N[ZTwoB[#1, xv1, vb1, vh1, xv2, vb2, vh2] & /@ ptsB, wp];
		predH = N[ZTwoH[#1, xv1, vb1, vh1, xv2, vb2, vh2] & /@ 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}
		]

(*   *)
printprogress[iteration_, chisq_, params_] := 
     Print[
    StringForm[
      "Iteration:`1` ChiSquared:`2`  {Position, ViabilityB, ViablityH, \
Position2, ViabilityB2, ViablityH2}:`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[]