function outputFile =generateSyntheticData(CONST)
%GENERATESYNTHETICDATA This function generates synthetic measurement data
%by creating and simulating a COMSOL model.

%  ========================================================================
fileName = [CONST.dataDirectory CONST.ComsolMeasurementTemplateFile];
disp(['   -> Loading COMSOL measurement model template from file "' fileName '"...']);
model = mphload(fileName, 'Measurement Model');
%  ========================================================================
disp('   -> Generating synthetic measurement points on virtual flight path...');
[measurementPoints, flightDetails] = generateFlightPath();
measurementSize = size(measurementPoints,1);
%  ========================================================================
if CONST.debugLevel >= 2
    disp('   -> Plotting virtual flight path...');
    figure('Name','Flight Path','NumberTitle','on');
    plotDomain(CONST.domainBds.xmin, CONST.domainBds.xmax, CONST.domainBds.ymin, CONST.domainBds.ymax, CONST.domainBds.zmin, CONST.domainBds.zmax);
    hold on; grid on;
    line(measurementPoints(:,1),measurementPoints(:,2),measurementPoints(:,3),'Color','black');
    scatter3(measurementPoints(:,1),measurementPoints(:,2),measurementPoints(:,3),15,'red','filled');
    title('Flight path measurement points');
end
%  ========================================================================
disp('   -> Writing measurement points to COMSOL model...');
for i=1:measurementSize
    mp = measurementPoints(i,:);
    dataId = ['measurement_pt' int2str(i)];
    model.geom('geom1').feature.create(dataId, 'Point');
    model.geom('geom1').feature(dataId).name(['Measurement Point ' int2str(i)]);
    model.geom('geom1').feature(dataId).set('p', {num2str(mp(1)); num2str(mp(2)); num2str(mp(3))});
end
model.geom('geom1').run();
%  ========================================================================
disp('   -> Writing measurement visualization elements to COMSOL model...');
model.result.create('pg_flight_path', 'PlotGroup3D');
model.result('pg_flight_path').name('Flight Path Measurements');
model.result('pg_flight_path').set('title', 'Measurement data on flight path');
for i=1 : flightDetails.numLanes
    dataId = ['measurement_pc' int2str(i)];
    model.result.dataset.create(dataId, 'ParCurve3D');
    model.result.dataset(dataId).name(['Flight Lane ' int2str(i)]);
    model.result.dataset(dataId).set('exprx', num2str(flightDetails.xStart+(i-1)*flightDetails.deltaX));
    model.result.dataset(dataId).set('expry', [num2str(flightDetails.yStart) ' + ' num2str(flightDetails.yEnd-flightDetails.yStart) '*s']);
    model.result.dataset(dataId).set('exprz', num2str(flightDetails.z));
    plotId = ['measurement_line' int2str(i)];
    model.result('pg_flight_path').feature.create(plotId, 'Line');
    model.result('pg_flight_path').feature(plotId).set('data', dataId);
    model.result('pg_flight_path').feature(plotId).set('linetype', 'tube');
    model.result('pg_flight_path').feature(plotId).set('radiusexpr', '10');
    model.result('pg_flight_path').feature(plotId).set('tuberadiusscaleactive', true);
    model.result('pg_flight_path').feature(plotId).set('tuberadiusscale', '10');
    model.result('pg_flight_path').feature(plotId).set('colortablesym', true);
    model.result('pg_flight_path').feature(plotId).set('colorlegend', 'off');
end
model.result.dataset.create('flight_surf', 'Surface');
model.result('pg_flight_path').feature.create('surf1', 'Surface');
model.result('pg_flight_path').feature.create('surf2', 'Surface');
model.result('pg_flight_path').feature('surf1').set('data', 'flight_surf');
model.result('pg_flight_path').feature('surf1').set('coloring', 'uniform');
model.result('pg_flight_path').feature('surf1').set('color', 'black');
model.result('pg_flight_path').feature('surf1').set('wireframe', true);
model.result('pg_flight_path').feature('surf2').set('data', 'flight_surf');
model.result('pg_flight_path').feature('surf2').set('coloring', 'uniform');
model.result('pg_flight_path').feature('surf2').set('color', 'gray');
%  ========================================================================
disp('   -> Computing mesh in COMSOL...');
model.mesh('mesh1').run();
%  ========================================================================
[pathstr, templateName, ext] = fileparts(CONST.ComsolMeasurementTemplateFile);
fileName = [CONST.outputDirectory 'PIPS-modified_' templateName '.mph'];
disp(['   -> Saving measurement model to file "' fileName '"...']);
model.sol('sol1').clearSolution;
mphsave(model, fileName);
%  ========================================================================
disp('   -> Simulating measurement model in COMSOL...');
model.sol('sol1').runAll;
u = mphgetu(model);
dof = length(u);
disp(['      -> degrees of freedom: ' int2str(dof)]);
%  ========================================================================
disp('   -> Extracting FEM data from COMSOL measurement model...');
% deactivate simulation for safety reasons (this is necessary):
model.sol('sol1').feature('s1').active(false); 
systemmatrix = mphmatrix(model,'sol1','out',{'K','L'});
meshinfo = mphxmeshinfo(model);
K = sparse(systemmatrix.K);
L = systemmatrix.L;
nodes = meshinfo.dofs.coords.';
%  ========================================================================
disp('   -> Assembling flight path look-up table and restriction operator...');
[T, measurementNodes] = assembleRestrictionOperator(measurementPoints, nodes, CONST.doubleComparisonTolerance);
%  ========================================================================
if CONST.debugLevel >=1
    disp('   -> Plotting restriction operator...');
    figure('Name','Restriction operator T for the measurement model','NumberTitle','on');
    spy(T,15);
    axis normal;
    title('Restriction operator T for the measurement model');
    xlabel('node index');
    ylabel('measurement index');
end
%  ========================================================================
disp('   -> Extracting measurement data on flight path...');
uGamma = T*u;
%  ========================================================================
disp('   -> Collecting and packing synthetic measurement data...');
syntheticData.measurementPoints = measurementPoints;
syntheticData.uGamma = uGamma;
syntheticData.COMSOL.K = K;
syntheticData.COMSOL.L = L;
syntheticData.COMSOL.u = u;
syntheticData.COMSOL.meshinfo = meshinfo;
syntheticData.T = T;
syntheticData.measurementNodes = measurementNodes;
%  ========================================================================
fileName = [CONST.outputDirectory 'PIPS-MeasData_from_' templateName '.mat'];
disp(['   -> Saving entire measurement data to file "' fileName '"...']);
save(fileName,'syntheticData');
%  ========================================================================
outputFile = fileName;
end