+ double complex amplitude = 0.0 + 0.0*I;
+ double complex lastamplitude = 0.0 + 0.0*I;
+
+ double complex probability = 1.0 + 0.0*I;
+
+ // the first measurement should be the identity
+ omega = 0;
+ for(i=0; i<N; i++) {
+ alpha[i] = 1; beta[i] = 0; gamma[i] = 0; delta[i] = 0;
+ }
+
+ for(j=0; j<numStabStates; j++) { // the kets
+ Gamma[j] *= measurepauli(N, &K[j], h[j], G[j], GBar[j], &Q[j], &D[j], &J[j], omega, gamma, beta);
+ }
+
+ double complex newamplitude;
+
+ #pragma omp parallel
+ {
+ #pragma omp for schedule(static) collapse(2)
+ for(i=0; i<numStabStates; i++) { // the bras
+ for(j=0; j<numStabStates; j++) {
+ newamplitude = innerproduct(N, K[j], h[j], G[j], GBar[j], Q[j], D[j], J[j], N, origK[i], origh[i], origG[i], origGBar[i], origQ[i], origD[i], origJ[i]);
+ amplitude = amplitude + conj(origGamma[i])*Gamma[j]*newamplitude;
+ }
+ }
+ }
+