% function out = model
%
% ITMPlusRing.m
%
% Model exported on Jun 24 2013, 16:24 by COMSOL 4.3.0.151.

u = zeros(1,2);
e = zeros(1,3);
mtime = 0;
deltat = 600;
totalt = 3600*24;
t = [0:deltat:totalt];
omega = 54/1000;
dnoverdt = (8.6)*10^(-6);
Ki = 0;
Kp = 0;
Kd = 0;
s_target = 1.5416e-05; %(m^-1) Corresponds to 5 W.

import com.comsol.model.*
import com.comsol.model.util.*

model = ModelUtil.create('Model');

model.modelPath('C:\Users\ligo\Desktop\SURF 2013');

model.name('ITMPlusRing.mph');

model.param.set('omega0', '54/1000');
model.param.set('r0', '170.065 / 1000');
model.param.set('th0', '199.61 / 1000');
model.param.set('I0', '109.2');
model.param.set('P0', '0');
model.param.set('Prh', '3');
model.param.set('deltat', num2str(deltat));
model.param.set('totalt', num2str(totalt));
model.param.set('z0', '4/100');
model.param.set('omegarh', '4/1000');
model.param.set('Irh', 'Prh/(2*pi*r0*0.00501326)');
model.param.set('mtime', num2str(mtime));

model.modelNode.create('mod1');

model.geom.create('geom1', 2);
model.geom('geom1').axisymmetric(true);
model.geom('geom1').feature.create('r1', 'Rectangle');
model.geom('geom1').feature.create('pt1', 'Point');
model.geom('geom1').feature.create('pt2', 'Point');
model.geom('geom1').feature.create('pt3', 'Point');
model.geom('geom1').feature.create('pc1', 'ParametricCurve');
model.geom('geom1').feature('r1').name('ITM');
model.geom('geom1').feature('r1').set('size', {'r0' 'th0'});
model.geom('geom1').feature('pt2').set('p', {'0'; '54 / 1000'});
model.geom('geom1').feature('pt3').set('p', {'54/1000'; '0'});
model.geom('geom1').feature('pc1').set('parmin', '-(6/1000)');
model.geom('geom1').feature('pc1').set('parmax', '(6/1000)');
model.geom('geom1').feature('pc1').set('coord', {'r0' 'th0 - (4/100) + s'});
model.geom('geom1').run;

% model.view.create('view2', 3);

model.material.create('mat1');
model.material('mat1').info.create('Composition');
model.material('mat1').propertyGroup('def').func.create('dL_solid_NIST_SRM_739_Type_I_1', 'Piecewise');
model.material('mat1').propertyGroup('def').func.create('CTE_solid_NIST_SRM_739_Type_I_1', 'Piecewise');
model.material('mat1').propertyGroup('def').func.create('k', 'Piecewise');
model.material('mat1').propertyGroup('def').func.create('res_solid_1', 'Piecewise');
model.material('mat1').propertyGroup('def').func.create('epsilon', 'Piecewise');
model.material('mat1').propertyGroup('def').func.create('alpha_solid_NIST_SRM_739_Type_I_1', 'Piecewise');
model.material('mat1').propertyGroup('def').func.create('C', 'Piecewise');
model.material('mat1').propertyGroup('def').func.create('HC', 'Piecewise');
model.material('mat1').propertyGroup('def').func.create('mu', 'Piecewise');
model.material('mat1').propertyGroup('def').func.create('sigma_solid_1', 'Piecewise');
model.material('mat1').propertyGroup('def').func.create('rho', 'Piecewise');
model.material('mat1').propertyGroup('def').func.create('TD', 'Piecewise');
model.material('mat1').propertyGroup('def').func.create('kappa', 'Piecewise');
model.material('mat1').propertyGroup.create('Enu', 'Young''s modulus and Poisson''s ratio');
model.material('mat1').propertyGroup('Enu').func.create('E', 'Piecewise');
model.material('mat1').propertyGroup('Enu').func.create('nu', 'Piecewise');

model.cpl.create('linproj1', 'LinearProjection', 'geom1');
model.cpl('linproj1').selection.set([1]);

model.physics.create('ht', 'HeatTransfer', 'geom1');
model.physics('ht').field('heatflux').field('J2');
model.physics('ht').feature.create('hf1', 'HeatFluxBoundary', 1);
model.physics('ht').feature('hf1').selection.set([2 5]);
model.physics('ht').feature.create('sar1', 'SurfaceToAmbientRadiation', 1);
model.physics('ht').feature('sar1').selection.all;
model.physics('ht').feature.create('hf2', 'HeatFluxBoundary', 1);
model.physics('ht').feature('hf2').selection.set([8]);

model.mesh.create('mesh1', 'geom1');
model.mesh('mesh1').feature.create('ftri1', 'FreeTri');

% model.result.table.create('evl3', 'Table');

% model.view('view1').axis.set('xmin', '-0.5940324068069458');
% model.view('view1').axis.set('xmax', '0.7936422824859619');
% model.view('view1').axis.set('ymin', '-0.4219888150691986');
% model.view('view1').axis.set('ymax', '0.7621188163757324');

model.material('mat1').name('Corning 7940 (fused silica) [solid,NIST SRM 739 - Type I]');
model.material('mat1').info('Composition').body('46.7 Si, 53.3 O (wt%)');
model.material('mat1').propertyGroup('def').func('dL_solid_NIST_SRM_739_Type_I_1').set('arg', 'T');
model.material('mat1').propertyGroup('def').func('dL_solid_NIST_SRM_739_Type_I_1').set('pieces', {'80.0' '1000.0' '9.433058E-5-1.712933E-6*T^1+7.374643E-9*T^2-1.113379E-11*T^3+8.078969E-15*T^4-2.330577E-18*T^5'});
model.material('mat1').propertyGroup('def').func('CTE_solid_NIST_SRM_739_Type_I_1').set('arg', 'T');
model.material('mat1').propertyGroup('def').func('CTE_solid_NIST_SRM_739_Type_I_1').set('pieces', {'80.0' '1000.0' '-1.712933E-6+1.474929E-8*T^1-3.340137E-11*T^2+3.231588E-14*T^3-1.165288E-17*T^4'});
model.material('mat1').propertyGroup('def').func('k').set('arg', 'T');
model.material('mat1').propertyGroup('def').func('k').set('pieces', {'5.0' '280.0' '0.09229359+0.003299521*T^1+4.343238E-5*T^2-2.258807E-7*T^3+3.049078E-10*T^4'; '280.0' '1400.0' '-0.9854696+0.01820907*T^1-5.290512E-5*T^2+7.552525E-8*T^3-5.008123E-11*T^4+1.311324E-14*T^5'; '1400.0' '1473.0' '-1798.509+4.192721*T^1-0.003251176*T^2+8.408144E-7*T^3'});
model.material('mat1').propertyGroup('def').func('res_solid_1').set('arg', 'T');
model.material('mat1').propertyGroup('def').func('res_solid_1').set('pieces', {'293.0' '823.0' 'exp(102.5761-0.4185428*T^1+8.035453E-4*T^2-7.491946E-7*T^3+2.719273E-10*T^4)'});
model.material('mat1').propertyGroup('def').func('epsilon').set('arg', 'T');
model.material('mat1').propertyGroup('def').func('epsilon').set('pieces', {'75.0' '312.0' '0.7243795+5.530512E-5*T^1'});
model.material('mat1').propertyGroup('def').func('alpha_solid_NIST_SRM_739_Type_I_1').set('arg', 'T');
model.material('mat1').propertyGroup('def').func('alpha_solid_NIST_SRM_739_Type_I_1').set('pieces', {'80.0' '1000.0' '-3.695793E-7+5.412844E-9*T^1-1.114924E-11*T^2+1.004725E-14*T^3-3.418476E-18*T^4'});
model.material('mat1').propertyGroup('def').func('C').set('arg', 'T');
model.material('mat1').propertyGroup('def').func('C').set('pieces', {'20.0' '115.0' '-41.6562+3.449789*T^1-0.002681788*T^2'; '115.0' '925.0' '61.49405+1.926794*T^1+0.004463629*T^2-1.704271E-5*T^3+1.88109E-8*T^4-7.062466E-12*T^5'; '925.0' '2000.0' '890.2414+0.3707459*T^1-1.107604E-4*T^2+3.139437E-8*T^3'});
model.material('mat1').propertyGroup('def').func('HC').set('arg', 'T');
model.material('mat1').propertyGroup('def').func('HC').set('pieces', {'20.0' '115.0' '-2.195319+0.2035589*T^1-1.778359E-4*T^2'; '115.0' '925.0' '3.394744+0.1154128*T^1+2.697214E-4*T^2-1.027028E-6*T^3+1.133063E-9*T^4-4.253438E-13*T^5'; '925.0' '2000.0' '53.46386+0.02165728*T^1-6.251109E-6*T^2+1.800145E-9*T^3'});
model.material('mat1').propertyGroup('def').func('mu').set('arg', 'T');
model.material('mat1').propertyGroup('def').func('mu').set('pieces', {'2.0' '170.0' '3.093791E10-2.513966E7*T^1+372325.4*T^2-2429.327*T^3+7.746996*T^4-0.009450154*T^5'; '170.0' '450.0' '3.087221E10-6676611.0*T^1+36600.0*T^2-35.83581*T^3'; '450.0' '1700.0' '2.914548E10+9388166.0*T^1-8173.552*T^2+3.504887*T^3-5.616985E-4*T^4'});
model.material('mat1').propertyGroup('def').func('sigma_solid_1').set('arg', 'T');
model.material('mat1').propertyGroup('def').func('sigma_solid_1').set('pieces', {'293.0' '823.0' '1/(exp(2.719273E-10*T^4-7.491946E-07*T^3+8.035453E-04*T^2-4.185428E-01*T+1.025761E+02))'});
model.material('mat1').propertyGroup('def').func('rho').set('arg', 'T');
model.material('mat1').propertyGroup('def').func('rho').set('pieces', {'80.0' '1000.0' '2219.39+0.01105999*T^1-4.708012E-5*T^2+6.980692E-8*T^3-4.995273E-11*T^4+1.430874E-14*T^5'});
model.material('mat1').propertyGroup('def').func('TD').set('arg', 'T');
model.material('mat1').propertyGroup('def').func('TD').set('pieces', {'82.0' '1039.0' '1.559138E-6-4.991949E-9*T^1+1.170891E-11*T^2-1.277613E-14*T^3+5.492523E-18*T^4'});
model.material('mat1').propertyGroup('def').func('kappa').set('arg', 'T');
model.material('mat1').propertyGroup('def').func('kappa').set('pieces', {'2.0' '170.0' '3.548804E10-7.249187E7*T^1+1081126.0*T^2-5664.153*T^3+10.51074*T^4'; '170.0' '950.0' '3.234657E10+2.830576E7*T^1-87346.28*T^2+175.3606*T^3-0.1614875*T^4+5.404807E-5*T^5'; '950.0' '1700.0' '1.341176E10+8.005859E7*T^1-86531.61*T^2+41.37606*T^3-0.007122425*T^4'});
model.material('mat1').propertyGroup('def').set('dL', 'dL_solid_NIST_SRM_739_Type_I_1(T[1/K])-dL_solid_NIST_SRM_739_Type_I_1(Tempref[1/K])');
model.material('mat1').propertyGroup('def').set('CTE', 'CTE_solid_NIST_SRM_739_Type_I_1(T[1/K])[1/K]');
model.material('mat1').propertyGroup('def').set('thermalconductivity', {'k(T[1/K])[W/(m*K)]' '0' '0' '0' 'k(T[1/K])[W/(m*K)]' '0' '0' '0' 'k(T[1/K])[W/(m*K)]'});
model.material('mat1').propertyGroup('def').set('resistivity', {'res_solid_1(T[1/K])[ohm*m]' '0' '0' '0' 'res_solid_1(T[1/K])[ohm*m]' '0' '0' '0' 'res_solid_1(T[1/K])[ohm*m]'});
model.material('mat1').propertyGroup('def').set('emissivity', 'epsilon(T[1/K])');
model.material('mat1').propertyGroup('def').set('thermalexpansioncoefficient', {'alpha_solid_NIST_SRM_739_Type_I_1(T[1/K])[1/K]+(Tempref-293[K])/(T-Tempref)*(alpha_solid_NIST_SRM_739_Type_I_1(T[1/K])[1/K]-alpha_solid_NIST_SRM_739_Type_I_1(Tempref[1/K])[1/K])' '0' '0' '0' 'alpha_solid_NIST_SRM_739_Type_I_1(T[1/K])[1/K]+(Tempref-293[K])/(T-Tempref)*(alpha_solid_NIST_SRM_739_Type_I_1(T[1/K])[1/K]-alpha_solid_NIST_SRM_739_Type_I_1(Tempref[1/K])[1/K])' '0' '0' '0' 'alpha_solid_NIST_SRM_739_Type_I_1(T[1/K])[1/K]+(Tempref-293[K])/(T-Tempref)*(alpha_solid_NIST_SRM_739_Type_I_1(T[1/K])[1/K]-alpha_solid_NIST_SRM_739_Type_I_1(Tempref[1/K])[1/K])'});
model.material('mat1').propertyGroup('def').set('heatcapacity', 'C(T[1/K])[J/(kg*K)]');
model.material('mat1').propertyGroup('def').set('HC', 'HC(T[1/K])[J/(mol*K)]');
model.material('mat1').propertyGroup('def').set('mu', 'mu(T[1/K])[Pa]');
model.material('mat1').propertyGroup('def').set('electricconductivity', {'sigma_solid_1(T[1/K])[S/m]' '0' '0' '0' 'sigma_solid_1(T[1/K])[S/m]' '0' '0' '0' 'sigma_solid_1(T[1/K])[S/m]'});
model.material('mat1').propertyGroup('def').set('density', 'rho(T[1/K])[kg/m^3]');
model.material('mat1').propertyGroup('def').set('TD', 'TD(T[1/K])[m^2/s]');
model.material('mat1').propertyGroup('def').set('kappa', 'kappa(T[1/K])[Pa]');
model.material('mat1').propertyGroup('def').addInput('temperature');
model.material('mat1').propertyGroup('def').addInput('strainreferencetemperature');
model.material('mat1').propertyGroup('Enu').func('E').set('arg', 'T');
model.material('mat1').propertyGroup('Enu').func('E').set('pieces', {'2.0' '170.0' '7.19219E10-7.955576E7*T^1+1200039.0*T^2-7440.215*T^3+21.72957*T^4-0.02394189*T^5'; '170.0' '1700.0' '6.906825E10+1.167539E7*T^1+11447.54*T^2-25.91953*T^3+0.01545627*T^4-3.022226E-6*T^5'});
model.material('mat1').propertyGroup('Enu').func('nu').set('arg', 'T');
model.material('mat1').propertyGroup('Enu').func('nu').set('pieces', {'2.0' '170.0' '0.1623859-3.46611E-4*T^1+5.42002E-6*T^2-2.78502E-8*T^3+4.878708E-11*T^4'; '170.0' '950.0' '0.1442531+2.285363E-4*T^1-9.253567E-7*T^2+1.832686E-9*T^3-1.635196E-12*T^4+5.369587E-16*T^5'; '950.0' '1700.0' '0.04194888+4.230553E-4*T^1-4.739473E-7*T^2+2.308844E-10*T^3-4.040605E-14*T^4'});
model.material('mat1').propertyGroup('Enu').set('youngsmodulus', 'E(T[1/K])[Pa]');
model.material('mat1').propertyGroup('Enu').set('poissonsratio', 'nu(T[1/K])');
model.material('mat1').propertyGroup('Enu').addInput('temperature');

model.cpl('linproj1').selection('srcvertex1').geom('geom1', 0);
model.cpl('linproj1').selection('srcvertex1').set([1]);
model.cpl('linproj1').selection('srcvertex2').geom('geom1', 0);
model.cpl('linproj1').selection('srcvertex2').set([5]);
model.cpl('linproj1').selection('srcvertex3').geom('geom1', 0);
model.cpl('linproj1').selection('srcvertex3').set([3]);
model.cpl('linproj1').selection('srcvertex4').geom('geom1', 0);
model.cpl('linproj1').selection('dstvertex1').geom('geom1', 0);
model.cpl('linproj1').selection('dstvertex1').set([1]);
model.cpl('linproj1').selection('dstvertex2').geom('geom1', 0);
model.cpl('linproj1').selection('dstvertex2').set([5]);
model.cpl('linproj1').selection('dstvertex3').geom('geom1', 0);

model.physics('ht').feature('hf1').set('q0', '(P0 / 0.00458044)*exp(-685.871*r^2)');
model.physics('ht').feature('hf1').set('qtot', '');
model.physics('ht').feature('hf1').name('Laser Heat Flux');
model.physics('ht').feature('sar1').set('epsilon_rad_mat', 'userdef');
model.physics('ht').feature('sar1').set('epsilon_rad', '0.9');
model.physics('ht').feature('hf2').set('HeatFluxType', 'TotalHeatFlux');
model.physics('ht').feature('hf2').set('q0', '');
model.physics('ht').feature('hf2').set('qtot', 'Prh');
model.physics('ht').feature('hf2').name('Ring Heater Flux');

model.mesh('mesh1').feature('size').set('hauto', 5);
model.mesh('mesh1').run;

% model.result.table('evl3').name('Evaluation 3D');
% model.result.table('evl3').comments('Interactive 3D values');

model.study.create('std1');
model.study('std1').feature.create('time', 'Transient');

model.sol.create('sol1');
model.sol('sol1').study('std1');
model.sol('sol1').attach('std1');
model.sol('sol1').feature.create('st1', 'StudyStep');
model.sol('sol1').feature.create('v1', 'Variables');
model.sol('sol1').feature.create('t1', 'Time');
model.sol('sol1').feature('t1').feature.create('fc1', 'FullyCoupled');
model.sol('sol1').feature('t1').feature.create('d1', 'Direct');
model.sol('sol1').feature('t1').feature.remove('fcDef');

% model.result.dataset.create('rev1', 'Revolve2D');
% model.result.create('pg1', 'PlotGroup3D');
% model.result('pg1').feature.create('surf1', 'Surface');
% model.result.create('pg2', 'PlotGroup2D');
% model.result('pg2').feature.create('con1', 'Contour');
% model.result('pg2').feature.create('arws1', 'ArrowSurface');
% model.result.create('pg3', 'PlotGroup1D');
% model.result('pg3').set('probetag', 'none');
% model.result('pg3').feature.create('ptgr1', 'PointGraph');
% model.result('pg3').feature('ptgr1').selection.set([4]);
% model.result('pg3').feature.create('ptgr2', 'PointGraph');
% model.result('pg3').feature('ptgr2').selection.set([1]);
% model.result.create('pg4', 'PlotGroup1D');
% model.result('pg4').set('probetag', 'none');
% model.result('pg4').feature.create('lngr1', 'LineGraph');
% model.result('pg4').feature('lngr1').selection.set([2 5]);
% model.result.export.create('plot1', 'Plot');
% model.result.export.create('plot2', 'Plot');
% model.result.export.create('plot3', 'Plot');
% model.result.export.create('plot4', 'Plot');

model.study('std1').feature('time').set('tlist', 'range(0,deltat,deltat)');

model.sol('sol1').attach('std1');
model.sol('sol1').feature('st1').name('Compile Equations: Time Dependent');
model.sol('sol1').feature('st1').set('studystep', 'time');
model.sol('sol1').feature('v1').set('control', 'time');
model.sol('sol1').feature('t1').set('control', 'time');
model.sol('sol1').feature('t1').set('tlist', 'range(0,deltat,deltat)');
model.sol('sol1').feature('t1').set('tstepsbdf', 'strict');
model.sol('sol1').feature('t1').set('initialstepbdf', 'deltat');
model.sol('sol1').feature('t1').set('initialstepbdfactive', true);
model.sol('sol1').feature('t1').set('maxstepbdf', 'deltat');
model.sol('sol1').feature('t1').set('maxstepbdfactive', true);
model.sol('sol1').feature('t1').set('maxorder', '2');
model.sol('sol1').feature('t1').feature('fc1').set('maxiter', '5');
model.sol('sol1').feature('t1').feature('fc1').set('jtech', 'once');
model.sol('sol1').feature('t1').feature('d1').set('linsolver', 'pardiso');
model.sol('sol1').runAll;

% model.result.dataset('rev1').set('startangle', '-90');
% model.result.dataset('rev1').set('revangle', '225');
% model.result('pg1').name('Temperature, 3D (ht)');
% model.result('pg1').feature('surf1').set('colortable', 'ThermalLight');
% model.result('pg2').name('Isothermal Contours (ht)');
% model.result('pg2').feature('con1').set('colortable', 'ThermalLight');
% model.result('pg2').feature('arws1').set('expr', {'ht.tfluxr' 'ht.tfluxz'});
% model.result('pg2').feature('arws1').set('descr', 'Total heat flux');
% model.result('pg2').feature('arws1').set('scale', '1.3378850992616805E-4');
% model.result('pg2').feature('arws1').set('color', 'gray');
% model.result('pg2').feature('arws1').set('scaleactive', false);
% model.result('pg3').set('xlabel', 'Time');
% model.result('pg3').set('ylabel', 'Temperature (K)');
% model.result('pg3').set('xlabelactive', false);
% model.result('pg3').set('ylabelactive', false);
% model.result('pg4').set('xlabel', 'r-coordinate (m)');
% model.result('pg4').set('ylabel', 'linproj1(T) (m*K)');
% model.result('pg4').set('xlabelactive', false);
% model.result('pg4').set('ylabelactive', false);
% model.result('pg4').feature('lngr1').set('expr', 'linproj1(T)');
% model.result('pg4').feature('lngr1').set('unit', 'm*K');
% model.result('pg4').feature('lngr1').set('descr', 'linproj1(T)');
% model.result('pg4').feature('lngr1').set('xdata', 'expr');
% model.result('pg4').feature('lngr1').set('xdataexpr', 'r');
% model.result('pg4').feature('lngr1').set('xdataunit', 'm');
% model.result('pg4').feature('lngr1').set('xdatadescr', 'r-coordinate');
% model.result.export('plot1').set('plotgroup', 'pg3');
% model.result.export('plot1').set('plot', 'ptgr1');
% model.result.export('plot1').set('filename', 'C:\Users\ligo\Desktop\SURF 2013\Testplot.txt');
% model.result.export('plot2').set('plotgroup', 'pg4');
% model.result.export('plot2').set('plot', 'lngr1');
% model.result.export('plot2').set('filename', 'C:\Users\ligo\Desktop\SURF 2013\Testplot 2.txt');
% model.result.export('plot3').set('plotgroup', 'pg4');
% model.result.export('plot3').set('plot', 'lngr1');
% model.result.export('plot3').set('filename', 'C:\Users\ligo\Desktop\SURF 2013\Gaussian Beam Data.txt');
% model.result.export('plot4').set('plotgroup', 'pg4');
% model.result.export('plot4').set('plot', 'lngr1');
% model.result.export('plot4').set('filename', 'C:\Users\ligo\Desktop\SURF 2013\PathIntegral.txt');

mtime = deltat;
r_0 = mphinterp(model,'linproj1(T)','coord',[0;0],'t',deltat);
r_54 = mphinterp(model,'linproj1(T)','coord',[54/1000;0],'t',deltat);
sdiff = (r_54 - r_0).*(dnoverdt/(2*omega^2));
s = sdiff;
e(1) = s - s_target;
u(1) = Ki*e(1)*deltat + Kp*e(1) + Kd*e(1)/deltat;
tic
for mtime = deltat:deltat:totalt
   model.param.set('mtime', num2str(mtime));
   model.param.set('actuation', num2str(u(1)));
   model.physics('ht').feature('hf2').set('qtot', 'Prh');
   model.sol('sol1').feature('v1').set('initmethod', 'sol');
   model.sol('sol1').feature('v1').set('initsol', 'sol1');
   model.sol('sol1').feature('v1').set('solnum', '2');
   model.sol('sol1').run;
   r_0(end+1) = mphinterp(model,'linproj1(T)','coord',[0;0],'t',deltat);
   r_54(end+1) = mphinterp(model,'linproj1(T)','coord',[54/1000;0],'t',deltat);
   sdiff(end+1) = (r_54(end) - r_0(end)).*(dnoverdt/(2*omega^2));
   s = sdiff(end);
   [u, e] = PID(s, u, e, Ki, Kp, Kd, deltat, s_target);
   toc
end
plot(t,sdiff)
grid on
xlabel('Time (s)')
ylabel('Defocus (m^{-1})')
title('Defocus at 1 beam radius as a function of time')