[code]
    do {
        photon_condition = ALIVE;
        acc_pathLength = 0;
        x = 0,y = 0,z = 0;
        pPacketCounter++;
        it_2 = 0;
   cosine_Zen_Cont = p1.getzenithDir();

        azimuth_cont = p1.getazimuthDir();

       cont_sintheta = sqrt(1.0 - (cosine_Zen_Cont*cosine_Zen_Cont));



        ux = cont_sintheta * cos(azimuth_cont);
        uy = cont_sintheta * sin(azimuth_cont);
        uz = cosine_Zen_Cont;



        do {


            randomNumber = p1.getRandNum();
            while ((randomNumber <= 0.0));
         varS_cont = p1.getdeltaS() / avgSkinScattCoeff;

            if (acc_pathLength + varS_cont >= photonPathLength[it_2]) {
                step_walk = photonPathLength[it_2] - acc_pathLength;


                x2 = x + (step_walk*ux);
                y2 = y + (step_walk*uy);
                z2 = z + (step_walk*uz);


                rhoPositionNow2 = p1.getRadialPos2(x2, y2, z2);

    /*SECTION FOR THE GIBBS'S SAMPLER WHERE WE COLLECT THE DISTANCES
    TRAVERSED FOR X, Y, AND Z WHILE THE PHOTON IS ALIVE; CALCULATE
    MEAN 1,2,3 -> SIGMA 1,2,3 -> COVARIANCE 1,2,3 -> MAYBE CORRELATION
                COEFFICIENT RHO BETWEEN ALL 3 VARIABLES*/



                //next we will create a radial index based off rho
                ir_2 = mc1.convertDoubToInt(rhoPositionNow2 / dr2);
if (ir_2 > NR) {ir_2 = NR;
                }

  Csph_pulsed[ir_2][it_2] += 1;
                it_2 += 1;
            }//if bracket
            x += ux * varS_cont;
            y += uy * varS_cont;
            z += uz * varS_cont;

            acc_pathLength += varS_cont;



            if (e == 0.0) {//isotropic scattering
                cosine_Zen_Cont = p1.getzenithDir();
            }
            else {
//get forward peak of scattering by using a random number between 0 and 1
    tempStore1 = (1.0 - (e*e)) / (1.0 - e + (2.0*e*phaseScattVal));
cosine_Zen_Cont = ((1.0 + e * e) - (tempStore1*tempStore1) / (2.0*e));
            }
cont_sintheta = sqrt(fabs(1.0 - (cosine_Zen_Cont*cosine_Zen_Cont)));


            cosine_azimuth = cos(azimuth_cont);
            if (azimuth_cont < PI) {
        sin_azimuth = sqrt(1.0 - (cosine_azimuth*cosine_azimuth));
            }//if
            else {
            sin_azimuth = -sqrt(1.0 - (cosine_azimuth*cosine_azimuth));

            }//else

            /*Make a new travel path for the photon via successive runs*/
            if (1.0 - fabs(uz) <= NEAR_PERPEN) {
                ux_New = cont_sintheta * cosine_azimuth;
                uy_New = cont_sintheta * sin_azimuth;
                uz_New = cosine_Zen_Cont * SIGN(uz);
            }//if

            else {
          tempStore1 = sqrt(1.0 - uz * uz);
          ux_New = cont_sintheta * (ux*uz*cosine_azimuth - uy *sin_azimuth)/(tempStore1 + (ux*cosine_Zen_Cont));

        uy_New = cont_sintheta * (uy*uz*cosine_azimuth + ux * sin_azimuth)/ (tempStore1 + (uy*cosine_Zen_Cont));

         uz_New = -1.0 * (cont_sintheta*cosine_azimuth*tempStore1) +  (uz*cosine_Zen_Cont);
            }
            ux = ux_New;
            uy = uy_New;
            uz = uz_New; //update trajectory





            pulsedXVal = (double*)calloc(photPacket, sizeof(double));
            pulsedYVal = (double*)calloc(photPacket, sizeof(double));
            pulsedZVal = (double*)calloc(photPacket, sizeof(double));
                up = 0;
             counterTotalX = photPacket + up;
             counterTotalY = photPacket + up;
             counterTotalZ = photPacket + up;
            for (int i = 0; i < counterTotalX; ++i) {

                        pulsedXVal[i] = x;
                        pulsedYVal[i] = y;
                        pulsedZVal[i] = z;

            }
            if (acc_pathLength >= maxPathLength) {

                photon_condition = DEAD;
            }

            ++up;
            counterTotalX += up;
            counterTotalY += up;
            counterTotalZ += up;
            pulsedXVal = (double*)realloc(pulsedXVal, counterTotalX);
            pulsedYVal = (double*)realloc(pulsedYVal, counterTotalY);
            pulsedZVal = (double*)realloc(pulsedZVal, counterTotalZ);


            //loopIncrementer++;



        } while (photon_condition == ALIVE);
        for (int i = 0; i < counterTotalX; ++i) {
            xCollect.push_back(pulsedXVal[i]);
        }
        for (int j = 0; j < counterTotalY; ++j) {
            yCollect.push_back(pulsedYVal[j]);
        }
        for (int k = 0; k < counterTotalZ; ++k) {
            zCollect.push_back(pulsedZVal[k]);
        }


        /*
        xCollectSum += x;
        yCollectSum += y;
        zCollectSum += z;
        */
    } while (pPacketCounter <= photPacket);

    //NEW
    for (it4 = xCollect.begin(); it4 != xCollect.end(); ++it4) {
        xCollectSum += *it4;

    }
    for (it5 = yCollect.begin(); it5 != yCollect.end(); ++it5) {
        yCollectSum += *it5;
        //yCollect.push_back(pulsedYVal[m]);
    }
    for (it6 = zCollect.begin(); it6 != zCollect.end(); ++it6) {
        //zCollect.push_back(pulsedZVal[n]);
        zCollectSum += *it6;
    }




    //NEW
    cout << "Pulsed simulation complete." << endl;
    averageX = xCollectSum / counterTotalX;
    averageY = yCollectSum / counterTotalY;
    averageZ = zCollectSum / counterTotalZ;

    xCollectArr = (double*)calloc(counterTotalX,sizeof(double));
    yCollectArr = (double*)calloc(counterTotalY, sizeof(double));
    zCollectArr= (double*)calloc(counterTotalZ, sizeof(double));
    for(int a; a<counterTotalX;++a){
        for (it4 = xCollect.begin(); it4 != xCollect.end(); ++it4) {
            xCollectArr[a] = *it4;
        }
    }
    for (int b; b<counterTotalY; ++b) {
        for (it5 = xCollect.begin(); it5 != xCollect.end(); ++it5) {
            yCollectArr[b] = *it5;
        }
    }
    for (int c; a<counterTotalZ; ++c) {
        for (it6 = xCollect.begin(); it6 != xCollect.end(); ++it6) {
            zCollectArr[c] = *it6;
        }
    }

    stdDevXRho = p1.getStdDev(averageX, counterTotalX, xCollectArr);
    stdDevYRho = p1.getStdDev(averageY, counterTotalY, yCollectArr);
    stdDevZRho = p1.getStdDev(averageZ, counterTotalZ, zCollectArr);

    xyCorrelCoeff = p1.getSimpleCorrelation(xCollectArr, yCollectArr, counterTotalX, averageX, averageY);
     xzCorrelCoeff = p1.getSimpleCorrelation(xCollectArr, zCollectArr, counterTotalX, averageX, averageZ);
    yzCorrelCoeff = p1.getSimpleCorrelation(yCollectArr, zCollectArr, counterTotalX, averageY, averageZ);
                                                                               //Above will be the correlation coefficients
   multCorrelR = p1.getMultipleCorrelation(xyCorrelCoeff, xzCorrelCoeff, yzCorrelCoeff);
    //multiple correlation coefficient

//goes with rest of declared variables
    double* chain_x1_1 = NULL;
    double* chain_y1_2 = NULL;
    double* chain_x2_1 = NULL;
    double* chain_z2_2 = NULL;
    double* chain_y3_1 = NULL;
    double* chain_z3_2 = NULL;

    chain_x1_1 = (double*)calloc(counterTotalX, sizeof(double));
    chain_y1_2 = (double*)calloc(counterTotalY, sizeof(double));
    chain_x2_1 = (double*)calloc(counterTotalX, sizeof(double));
    chain_z2_2 = (double*)calloc(counterTotalZ, sizeof(double));
    chain_y3_1 = (double*)calloc(counterTotalY, sizeof(double));
    chain_z3_2 = (double*)calloc(counterTotalZ, sizeof(double));


          if (chain_x1_1 == NULL || chain_y1_2 == NULL || chain_x2_1 == NULL || chain_z2_2 == NULL || chain_y3_1 == NULL || chain_z3_2 == NULL) {
        cerr << "Memory allocation failed(2)" << endl;
    }


 for (int i = 0; i<counterTotalX; i++) {
        chain_x1_1[0] = 0;
        chain_y1_2[0] = 0;

        chain_x2_1[0] = 0;
        chain_z2_2[0] = 0;

        chain_y3_1[0] = 0;
        chain_z3_2[0] = 0;


        chain_x1_1[i + 1] = p1.getGibbsConditionalSample(chain_y1_2[i], averageX, averageY, stdDevXRho, stdDevYRho, xyCorrelCoeff);
         chain_y1_2[i + 1] = p1.getGibbsConditionalSample(chain_x1_1[i + 1], averageY, averageX, stdDevYRho, stdDevXRho, xyCorrelCoeff);

         [/code]