double Iq(double q, double case_num,
    double N[], double Phi[], double v[], double L[], double b[],
    double Kab, double Kac, double Kad,
    double Kbc, double Kbd, double Kcd
    );

double Iq(double q, double case_num,
    double N[], double Phi[], double v[], double L[], double b[],
    double Kab, double Kac, double Kad,
    double Kbc, double Kbd, double Kcd
    ) {
  int icase = (int)case_num;

#if 0  // Sasview defaults
  if (icase <= 1) {
    N[0]=N[1]=1000.0;
    Phi[0]=Phi[1]=0.0000001;
    Kab=Kac=Kad=Kbc=Kbd=-0.0004;
    L[0]=L[1]=1e-12;
    v[0]=v[1]=100.0;
    b[0]=b[1]=5.0;
  } else if (icase <= 4) {
    Phi[0]=0.0000001;
    Kab=Kac=Kad=-0.0004;
    L[0]=1e-12;
    v[0]=100.0;
    b[0]=5.0;
  }
#else
  if (icase <= 1) {
    N[0]=N[1]=0.0;
    Phi[0]=Phi[1]=0.0;
    Kab=Kac=Kad=Kbc=Kbd=0.0;
    L[0]=L[1]=L[3];
    v[0]=v[1]=v[3];
    b[0]=b[1]=0.0;
  } else if (icase <= 4) {
    N[0] = 0.0;
    Phi[0]=0.0;
    Kab=Kac=Kad=0.0;
    L[0]=L[3];
    v[0]=v[3];
    b[0]=0.0;
  }
#endif

  const double Xa = q*q*b[0]*b[0]*N[0]/6.0;
  const double Xb = q*q*b[1]*b[1]*N[1]/6.0;
  const double Xc = q*q*b[2]*b[2]*N[2]/6.0;
  const double Xd = q*q*b[3]*b[3]*N[3]/6.0;

  // limit as Xa goes to 0 is 1
  const double Pa = Xa==0 ? 1.0 : -expm1(-Xa)/Xa;
  const double Pb = Xb==0 ? 1.0 : -expm1(-Xb)/Xb;
  const double Pc = Xc==0 ? 1.0 : -expm1(-Xc)/Xc;
  const double Pd = Xd==0 ? 1.0 : -expm1(-Xd)/Xd;

  // limit as Xa goes to 0 is 1
  const double Paa = Xa==0 ? 1.0 : 2.0*(1.0-Pa)/Xa;
  const double Pbb = Xb==0 ? 1.0 : 2.0*(1.0-Pb)/Xb;
  const double Pcc = Xc==0 ? 1.0 : 2.0*(1.0-Pc)/Xc;
  const double Pdd = Xd==0 ? 1.0 : 2.0*(1.0-Pd)/Xd;


  // Note: S0ij only defined for copolymers; otherwise set to zero
  // 0: C/D     binary mixture
  // 1: C-D     diblock copolymer
  // 2: B/C/D   ternery mixture
  // 3: B/C-D   binary mixture,1 homopolymer, 1 diblock copolymer
  // 4: B-C-D   triblock copolymer
  // 5: A/B/C/D quaternary mixture
  // 6: A/B/C-D ternery mixture, 2 homopolymer, 1 diblock copolymer
  // 7: A/B-C-D binary mixture, 1 homopolymer, 1 triblock copolymer
  // 8: A-B/C-D binary mixture, 2 diblock copolymer
  // 9: A-B-C-D tetra-block copolymer
#if 0
  const double S0aa = icase<5
                      ? 1.0 : N[0]*Phi[0]*v[0]*Paa;
  const double S0bb = icase<2
                      ? 1.0 : N[1]*Phi[1]*v[1]*Pbb;
  const double S0cc = N[2]*Phi[2]*v[2]*Pcc;
  const double S0dd = N[3]*Phi[3]*v[3]*Pdd;
  const double S0ab = icase<8
                      ? 0.0 : sqrt(N[0]*v[0]*Phi[0]*N[1]*v[1]*Phi[1])*Pa*Pb;
  const double S0ac = icase<9
                      ? 0.0 : sqrt(N[0]*v[0]*Phi[0]*N[2]*v[2]*Phi[2])*Pa*Pc*exp(-Xb);
  const double S0ad = icase<9
                      ? 0.0 : sqrt(N[0]*v[0]*Phi[0]*N[3]*v[3]*Phi[3])*Pa*Pd*exp(-Xb-Xc);
  const double S0bc = (icase!=4 && icase!=7 && icase!= 9)
                      ? 0.0 : sqrt(N[1]*v[1]*Phi[1]*N[2]*v[2]*Phi[2])*Pb*Pc;
  const double S0bd = (icase!=4 && icase!=7 && icase!= 9)
                      ? 0.0 : sqrt(N[1]*v[1]*Phi[1]*N[3]*v[3]*Phi[3])*Pb*Pd*exp(-Xc);
  const double S0cd = (icase==0 || icase==2 || icase==5)
                      ? 0.0 : sqrt(N[2]*v[2]*Phi[2]*N[3]*v[3]*Phi[3])*Pc*Pd;
#else  // sasview equivalent
//printf("Xc=%g, S0cc=%g*%g*%g*%g\n",Xc,N[2],Phi[2],v[2],Pcc);
  double S0aa = N[0]*Phi[0]*v[0]*Paa;
  double S0bb = N[1]*Phi[1]*v[1]*Pbb;
  double S0cc = N[2]*Phi[2]*v[2]*Pcc;
  double S0dd = N[3]*Phi[3]*v[3]*Pdd;
  double S0ab = sqrt(N[0]*v[0]*Phi[0]*N[1]*v[1]*Phi[1])*Pa*Pb;
  double S0ac = sqrt(N[0]*v[0]*Phi[0]*N[2]*v[2]*Phi[2])*Pa*Pc*exp(-Xb);
  double S0ad = sqrt(N[0]*v[0]*Phi[0]*N[3]*v[3]*Phi[3])*Pa*Pd*exp(-Xb-Xc);
  double S0bc = sqrt(N[1]*v[1]*Phi[1]*N[2]*v[2]*Phi[2])*Pb*Pc;
  double S0bd = sqrt(N[1]*v[1]*Phi[1]*N[3]*v[3]*Phi[3])*Pb*Pd*exp(-Xc);
  double S0cd = sqrt(N[2]*v[2]*Phi[2]*N[3]*v[3]*Phi[3])*Pc*Pd;
switch(icase){
  case 0:
    S0aa=0.000001;
    S0ab=0.000002;
    S0ac=0.000003;
    S0ad=0.000004;
    S0bb=0.000005;
    S0bc=0.000006;
    S0bd=0.000007;
    S0cd=0.000008;
    break;
  case 1:
    S0aa=0.000001;
    S0ab=0.000002;
    S0ac=0.000003;
    S0ad=0.000004;
    S0bb=0.000005;
    S0bc=0.000006;
    S0bd=0.000007;
    break;
  case 2:
    S0aa=0.000001;
    S0ab=0.000002;
    S0ac=0.000003;
    S0ad=0.000004;
    S0bc=0.000005;
    S0bd=0.000006;
    S0cd=0.000007;
    break;
  case 3:
    S0aa=0.000001;
    S0ab=0.000002;
    S0ac=0.000003;
    S0ad=0.000004;
    S0bc=0.000005;
    S0bd=0.000006;
    break;
  case 4:
    S0aa=0.000001;
    S0ab=0.000002;
    S0ac=0.000003;
    S0ad=0.000004;
    break;
  case 5:
    S0ab=0.000001;
    S0ac=0.000002;
    S0ad=0.000003;
    S0bc=0.000004;
    S0bd=0.000005;
    S0cd=0.000006;
    break;
  case 6:
    S0ab=0.000001;
    S0ac=0.000002;
    S0ad=0.000003;
    S0bc=0.000004;
    S0bd=0.000005;
    break;
  case 7:
    S0ab=0.000001;
    S0ac=0.000002;
    S0ad=0.000003;
    break;
  case 8:
    S0ac=0.000001;
    S0ad=0.000002;
    S0bc=0.000003;
    S0bd=0.000004;
    break;
  default : //case 9:
    break;
  }
#endif

  // eq 12a: \kappa_{ij}^F = \chi_{ij}^F - \chi_{i0}^F - \chi_{j0}^F
  const double Kaa = 0.0;
  const double Kbb = 0.0;
  const double Kcc = 0.0;
  //const double Kdd = 0.0;
  const double Zaa = Kaa - Kad - Kad;
  const double Zab = Kab - Kad - Kbd;
  const double Zac = Kac - Kad - Kcd;
  const double Zbb = Kbb - Kbd - Kbd;
  const double Zbc = Kbc - Kbd - Kcd;
  const double Zcc = Kcc - Kcd - Kcd;
//printf("Za: %10.5g %10.5g %10.5g\n", Zaa, Zab, Zac);
//printf("Zb: %10.5g %10.5g %10.5g\n", Zab, Zbb, Zbc);
//printf("Zc: %10.5g %10.5g %10.5g\n", Zac, Zbc, Zcc);

  // T = inv(S0)
  const double DenT = (- S0ac*S0bb*S0ac + S0ab*S0bc*S0ac + S0ac*S0ab*S0bc
                       - S0aa*S0bc*S0bc - S0ab*S0ab*S0cc + S0aa*S0bb*S0cc);
  const double T11 = (-S0bc*S0bc + S0bb*S0cc)/DenT;
  const double T12 = ( S0ac*S0bc - S0ab*S0cc)/DenT;
  const double T13 = (-S0ac*S0bb + S0ab*S0bc)/DenT;
  const double T22 = (-S0ac*S0ac + S0aa*S0cc)/DenT;
  const double T23 = ( S0ac*S0ab - S0aa*S0bc)/DenT;
  const double T33 = (-S0ab*S0ab + S0aa*S0bb)/DenT;

//printf("T1: %10.5g %10.5g %10.5g\n", T11, T12, T13);
//printf("T2: %10.5g %10.5g %10.5g\n", T12, T22, T23);
//printf("T3: %10.5g %10.5g %10.5g\n", T13, T23, T33);

  // eq 18e: m = 1/(S0_{dd} - s0^T inv(S0) s0)
  const double ZZ = S0ad*(T11*S0ad + T12*S0bd + T13*S0cd)
                  + S0bd*(T12*S0ad + T22*S0bd + T23*S0cd)
                  + S0cd*(T13*S0ad + T23*S0bd + T33*S0cd);

  const double m=1.0/(S0dd-ZZ);

  // eq 18d: Y = inv(S0)s0 + e
  const double Y1 = T11*S0ad + T12*S0bd + T13*S0cd + 1.0;
  const double Y2 = T12*S0ad + T22*S0bd + T23*S0cd + 1.0;
  const double Y3 = T13*S0ad + T23*S0bd + T33*S0cd + 1.0;

  // N = mYY^T + \kappa^F
  const double N11 = m*Y1*Y1 + Zaa;
  const double N12 = m*Y1*Y2 + Zab;
  const double N13 = m*Y1*Y3 + Zac;
  const double N22 = m*Y2*Y2 + Zbb;
  const double N23 = m*Y2*Y3 + Zbc;
  const double N33 = m*Y3*Y3 + Zcc;

//printf("N1: %10.5g %10.5g %10.5g\n", N11, N12, N13);
//printf("N2: %10.5g %10.5g %10.5g\n", N12, N22, N23);
//printf("N3: %10.5g %10.5g %10.5g\n", N13, N23, N33);
//printf("S0a: %10.5g %10.5g %10.5g\n", S0aa, S0ab, S0ac);
//printf("S0b: %10.5g %10.5g %10.5g\n", S0ab, S0bb, S0bc);
//printf("S0c: %10.5g %10.5g %10.5g\n", S0ac, S0bc, S0cc);

  // M = I + S0 N
  const double Maa = N11*S0aa + N12*S0ab + N13*S0ac + 1.0;
  const double Mab = N11*S0ab + N12*S0bb + N13*S0bc;
  const double Mac = N11*S0ac + N12*S0bc + N13*S0cc;
  const double Mba = N12*S0aa + N22*S0ab + N23*S0ac;
  const double Mbb = N12*S0ab + N22*S0bb + N23*S0bc + 1.0;
  const double Mbc = N12*S0ac + N22*S0bc + N23*S0cc;
  const double Mca = N13*S0aa + N23*S0ab + N33*S0ac;
  const double Mcb = N13*S0ab + N23*S0bb + N33*S0bc;
  const double Mcc = N13*S0ac + N23*S0bc + N33*S0cc + 1.0;
//printf("M1: %10.5g %10.5g %10.5g\n", Maa, Mab, Mac);
//printf("M2: %10.5g %10.5g %10.5g\n", Mba, Mbb, Mbc);
//printf("M3: %10.5g %10.5g %10.5g\n", Mca, Mcb, Mcc);

  // Q = inv(M) = inv(I + S0 N)
  const double DenQ = (+ Maa*Mbb*Mcc - Maa*Mbc*Mcb - Mab*Mba*Mcc
                       + Mab*Mbc*Mca + Mac*Mba*Mcb - Mac*Mbb*Mca);

  const double Q11 = ( Mbb*Mcc - Mbc*Mcb)/DenQ;
  const double Q12 = (-Mab*Mcc + Mac*Mcb)/DenQ;
  const double Q13 = ( Mab*Mbc - Mac*Mbb)/DenQ;
  //const double Q21 = (-Mba*Mcc + Mbc*Mca)/DenQ;
  const double Q22 = ( Maa*Mcc - Mac*Mca)/DenQ;
  const double Q23 = (-Maa*Mbc + Mac*Mba)/DenQ;
  //const double Q31 = ( Mba*Mcb - Mbb*Mca)/DenQ;
  //const double Q32 = (-Maa*Mcb + Mab*Mca)/DenQ;
  const double Q33 = ( Maa*Mbb - Mab*Mba)/DenQ;

//printf("Q1: %10.5g %10.5g %10.5g\n", Q11, Q12, Q13);
//printf("Q2: %10.5g %10.5g %10.5g\n", Q21, Q22, Q23);
//printf("Q3: %10.5g %10.5g %10.5g\n", Q31, Q32, Q33);
  // eq 18c: inv(S) = inv(S0) + mYY^T + \kappa^F
  // eq A1 in the appendix
  // To solve for S, use:
  //      S = inv(inv(S^0) + N) inv(S^0) S^0
  //        = inv(S^0 inv(S^0) + N) S^0
  //        = inv(I + S^0 N) S^0
  //        = Q S^0
  const double S11 = Q11*S0aa + Q12*S0ab + Q13*S0ac;
  const double S12 = Q12*S0aa + Q22*S0ab + Q23*S0ac;
  const double S13 = Q13*S0aa + Q23*S0ab + Q33*S0ac;
  const double S22 = Q12*S0ab + Q22*S0bb + Q23*S0bc;
  const double S23 = Q13*S0ab + Q23*S0bb + Q33*S0bc;
  const double S33 = Q13*S0ac + Q23*S0bc + Q33*S0cc;
  // If the full S is needed...it isn't since Ldd = (rho_d - rho_d) = 0 below
  //const double S14=-S11-S12-S13;
  //const double S24=-S12-S22-S23;
  //const double S34=-S13-S23-S33;
  //const double S44=S11+S22+S33 + 2.0*(S12+S13+S23);

  // eq 12 of Akcasu, 1990: I(q) = L^T S L
  // Note: eliminate cases without A and B polymers by setting Lij to 0
  // Note: 1e-13 to convert from fm to cm for scattering length
  const double sqrt_Nav=sqrt(6.022045e+23) * 1.0e-13;
  const double Lad = icase<5 ? 0.0 : (L[0]/v[0] - L[3]/v[3])*sqrt_Nav;
  const double Lbd = icase<2 ? 0.0 : (L[1]/v[1] - L[3]/v[3])*sqrt_Nav;
  const double Lcd = (L[2]/v[2] - L[3]/v[3])*sqrt_Nav;

  const double result=Lad*Lad*S11 + Lbd*Lbd*S22 + Lcd*Lcd*S33
                    + 2.0*(Lad*Lbd*S12 + Lbd*Lcd*S23 + Lad*Lcd*S13);

  return result;

}