Index: issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/README.md
===================================================================
--- issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/README.md	(revision 26975)
+++ issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/README.md	(revision 26975)
@@ -0,0 +1,15 @@
+# Function List
+
+* analyzeCalvingFront.m				--		Interpolate the transient solutions to the calving front positions of the given flowlines
+* averageOverTime.m					--		For a given range, average transient data in time
+* computeGrad.m						--		Compute the gradient of a given field
+* extractTransientSolutions.m		--		Load a model with transient solutions, and put each of the items into individual array
+* interpZeroPos.m						--		Find zero positions (y=0) of the given curve defined by (X,Y)
+* integrateOverDomain.m				--		integrate a given variable over the whole domain, can set a mask
+* integrateOverTime.m				--		integrate a time series over the given time steps
+* inverseHistograms.m				--		Compute the inverse quantile function
+* npsd.m									--		Compute normalized power spectral density
+* projectToFlowlines.m				--		Project data from a mesh to a given flowline
+* psd.m									--		Compute power spectral density 
+* rescalegradientNan.m				--		Replace Nan by 0 in the cost function, then rescale it
+* wfDistance.m							--		Compute Wasserstein-Fourier distance
Index: issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/analyzeCalvingFront.m
===================================================================
--- issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/analyzeCalvingFront.m	(revision 26975)
+++ issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/analyzeCalvingFront.m	(revision 26975)
@@ -0,0 +1,62 @@
+function [flowline, icemaskFL] = analyzeCalvingFront(md, flowline, transientSolutions)
+%analyzeCalvingFront - to find calving front position at the given flowline,
+%                   if the given model contains transient solutions, then
+%                   return a time dependent calving front with
+%                   [Xmain, positionx, positiony, time], icemask, calving rate,
+%                   melting rate and velocity magnitude, sigma, thickness
+% Author: Cheng Gong
+% Last modified: 2020-09-25
+
+transient = isfield(md.results,'TransientSolution');
+
+if transient
+    % extract data from model
+    icemaskFL = InterpFromMeshToMesh2d(md.mesh.elements,md.mesh.x,md.mesh.y,transientSolutions.ice_levelset,flowline.x,flowline.y);
+    
+    % solve for the calving front coordinates
+    cf = interpZeroPos([flowline.Xmain(:), flowline.x, flowline.y], icemaskFL);
+    positionx = cf(:, 2);
+    positiony = cf(:, 3);
+    
+    % Calving rate C
+    cRate = InterpFromMeshToMesh2d(md.mesh.elements,md.mesh.x,md.mesh.y,...
+		 transientSolutions.calvingRate,positionx, positiony);
+    cRate = diag(cRate); % not very efficient, but works
+    
+    % melting rate M
+    mRateTemp = InterpFromMeshToMesh2d(md.mesh.elements,md.mesh.x,md.mesh.y,...
+		 transientSolutions.meltingRate,positionx, positiony);
+    mRate = diag(mRateTemp); % not very efficient, but works
+    
+    % velocity
+    vel = InterpFromMeshToMesh2d(md.mesh.elements,md.mesh.x,md.mesh.y,...
+		 transientSolutions.vel, positionx, positiony);
+    vel = diag(vel);
+    
+    % sigmaVM
+    SigmaVM = InterpFromMeshToMesh2d(md.mesh.elements,md.mesh.x,md.mesh.y,...
+		 transientSolutions.SigmaVM, positionx, positiony);
+    SigmaVM = diag(SigmaVM);
+    
+    % calving front
+    calvingFront = [cf, transientSolutions.time(:)];
+    
+    % solutions along the flowline
+    SigmaVM_FL = InterpFromMeshToMesh2d(md.mesh.elements,md.mesh.x,md.mesh.y,...
+		 transientSolutions.SigmaVM,flowline.x, flowline.y);
+    
+    % Thickness along the flowlines
+    thickness = InterpFromMeshToMesh2d(md.mesh.elements,md.mesh.x,md.mesh.y,...
+		 transientSolutions.thickness, flowline.x, flowline.y);
+    
+    % save to flowline
+    flowline.calvingFront = calvingFront;
+    flowline.cRate = cRate;
+    flowline.mRate = mRate;
+    flowline.vel = vel;
+    flowline.SigmaVM = SigmaVM;
+    flowline.SigmaVM_FL = SigmaVM_FL;
+    flowline.thickness = thickness;
+else
+    error('Calving front detectiong for the steady state is not implemented yet!')
+end
Index: issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/averageOverTime.m
===================================================================
--- issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/averageOverTime.m	(revision 26975)
+++ issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/averageOverTime.m	(revision 26975)
@@ -0,0 +1,86 @@
+function averagedData = averageOverTime(data, time, startP, endP)
+%averageOverTime - compute the time averaged value of data in the range [startP, endP]
+%
+%   data: time dependent data,
+%   time: has same number of columns as data, sorted in ascending order
+%   startP: start time point for averaging, can be at any time point, or
+%   even not coincide with the given points in 'time'
+%   endP: end time point for averaging. if this is the same as the startP,
+%   or not given, then the output will be the linear interpolation of the
+%   data at startP.
+%
+% if out of the given time range, use constant extrapolation,
+% if within the range, do linear interpolation.
+%
+%   Author: Cheng Gong
+%   Last modified: 2020-09-09
+if nargin < 4
+    endP = startP;
+end
+
+Nt = length(time);
+if Nt ~= size(data,2)
+    error('The data and time does not match!');
+end
+
+if startP > endP
+    error('The given start point has to be smaller than the end point!');
+elseif startP == endP
+    % extrapolation if earlier than the first time point
+    if startP <= time(1)
+        i1 = 1;
+        i2 = 2;
+    % extrapolation if later than the last time point
+    elseif startP >= time(Nt)
+        i1 = Nt-1;
+        i2 = Nt;
+    else
+        pos = (time > startP);
+        i2 = find(pos, 1, 'first');
+        i1 = i2 - 1;
+    end
+    % linear interpolation 
+    averagedData = data(:, i1) + (data(:, i1)-data(:, i2)).*(startP-time(i1))./(time(i1)-time(i2));
+else
+    % find the first and last ID in the time series within the given range
+    pos = ((time>=startP) & (time<=endP));
+    firstId = find(pos, 1, 'first');
+    lastId = find(pos, 1, 'last');
+    if isempty(firstId) | isempty(lastId)
+		 averagedData = nan(size(data,1), 1);
+       % error('Time points out of range');
+	 else
+		 % computed the integral with trapzoidal rule
+   	 integData = zeros(size(data,1), 1);
+   	 for i = firstId:lastId-1
+   	     integData = integData + (data(:,i+1) + data(:,i)) .* (time(i+1) - time(i)).*0.5;
+   	 end
+   	 
+   	 % special treatment for the first and last time step, which are not
+   	 % complete steps.
+   	 if firstId == 1
+   	     % extrapolation
+   	     integData = integData + data(:,1) .* (time(1)-startP);
+   	 else
+   	     % interpolation
+   	     integData = integData + 0.5.*(time(firstId)-startP)./ ...
+   	         (time(firstId)-time(firstId-1)) .* ...
+   	         ((time(firstId)-startP) .* data(:,firstId-1) + ...
+   	         (time(firstId)+startP-2*time(firstId-1)).*data(:,firstId));
+   	 end
+   	 
+   	 if lastId == Nt
+   	     % extrapolation
+   	     integData = integData + data(:,Nt) .* (endP-time(Nt));
+   	 else
+   	     % interpolation
+   	     integData = integData + 0.5.*(time(lastId)-endP)./ ...
+   	         (time(lastId+1)-time(lastId)) .* ...
+   	         ((time(lastId)-endP) .* data(:,lastId+1) + ...
+   	         (time(lastId)+endP-2*time(lastId+1)).*data(:,lastId));
+   	 end
+   	 
+   	 averagedData = integData ./ (endP-startP);
+	end
+end
+
Index: issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/computeGrad.m
===================================================================
--- issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/computeGrad.m	(revision 26975)
+++ issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/computeGrad.m	(revision 26975)
@@ -0,0 +1,45 @@
+function [gradx, grady]=computeGrad(index,x,y,field)
+%COMPUTEHESSIAN - compute the gradient from a field
+
+%some variables
+numberofnodes=length(x);
+numberofelements=size(index,1);
+numberoftime = size(field, 2);
+
+%some checks
+if (length(field)~=numberofnodes) && (length(field)~=numberofelements)
+    error('ComputeHessian error message: the given field size not supported yet');
+end
+%initialization
+line=index(:);
+linesize=3*numberofelements;
+
+%get areas and nodal functions coefficients N(x,y)=alpha x + beta y + gamma
+[alpha, beta]=GetNodalFunctionsCoeff(index,x,y);
+areas=GetAreas(index,x,y);
+
+%compute weights that hold the volume of all the element holding the node i
+weights=sparse(line,ones(linesize,1),repmat(areas,3,1),numberofnodes,1);
+
+%compute field on nodes if on elements
+if length(field)==numberofelements
+    field=sparse(line,ones(linesize,1),repmat(areas.*field,3,1),numberofnodes,1)./weights ;
+end
+
+%Compute gradient for each element
+if numberoftime == 1
+    grad_elx=sum(field(index).*alpha,2);
+    grad_ely=sum(field(index).*beta,2);
+else
+    grad_elx = zeros(numberofelements,numberoftime);
+    grad_ely = zeros(numberofelements,numberoftime);    
+    for i = 1:3
+        grad_elx = grad_elx + field(index(:,i), :).*alpha(:,i);
+        grad_ely = grad_ely + field(index(:,i), :).*beta(:,i);
+    end
+end
+%Compute gradient for each node (average of the elements around)
+gradx=sparse(repmat(line,1,numberoftime),cumsum(ones(linesize,numberoftime),2),repmat(areas.*grad_elx,3,1),numberofnodes,numberoftime);
+grady=sparse(repmat(line,1,numberoftime),cumsum(ones(linesize,numberoftime),2),repmat(areas.*grad_ely,3,1),numberofnodes,numberoftime);
+gradx=gradx./weights;
+grady=grady./weights;
Index: issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/extractTransientSolutions.m
===================================================================
--- issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/extractTransientSolutions.m	(revision 26975)
+++ issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/extractTransientSolutions.m	(revision 26975)
@@ -0,0 +1,14 @@
+function [transientSolutions]=extractTransientSolutions(md)
+%extractTransientSolutions - take the transient solutions out from model
+%                            and put each of them into an individual array.
+transientSolutions.time = cell2mat({md.results.TransientSolution(:).time});
+transientSolutions.vx = cell2mat({md.results.TransientSolution(:).Vx});
+transientSolutions.vy = cell2mat({md.results.TransientSolution(:).Vy});
+transientSolutions.vel = cell2mat({md.results.TransientSolution(:).Vel});
+transientSolutions.volume = cell2mat({md.results.TransientSolution(:).IceVolume});
+transientSolutions.thickness = cell2mat({md.results.TransientSolution(:).Thickness});
+transientSolutions.SigmaVM = cell2mat({md.results.TransientSolution(:).SigmaVM});
+transientSolutions.smb = cell2mat({md.results.TransientSolution(:).SmbMassBalance});
+transientSolutions.ice_levelset = cell2mat({md.results.TransientSolution(:).MaskIceLevelset});
+transientSolutions.calvingRate = cell2mat({md.results.TransientSolution(:).CalvingCalvingrate});
+transientSolutions.meltingRate = cell2mat({md.results.TransientSolution(:).CalvingMeltingrate});
Index: issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/fillInNan.m
===================================================================
--- issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/fillInNan.m	(revision 26975)
+++ issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/fillInNan.m	(revision 26975)
@@ -0,0 +1,26 @@
+function newdata = fillInNan(md, data, varargin)
+%fillInNan - use the mean of surrouding data to iteratidatay fill in the Nan in data
+%		data should have the same size as md.mesh.x
+%
+% Author: Cheng Gong
+% Last modified: 2021-12-08
+nanvflag = find(isnan(data));
+NNanv = length(nanvflag);
+
+options    = pairoptions(varargin{:});
+maxiter    = getfieldvalue(options,'maxiter', 10);
+count = 1;
+
+while((NNanv>0) & (count<=maxiter))
+   disp(['Iter ', num2str(count), ': found ', num2str(NNanv), ' Nan in the initial data, fill in them with the mean of their surroudings.']);
+   for i = 1:NNanv
+      [eleid, ~] = find(md.mesh.elements == nanvflag(i));
+      nodeid = md.mesh.elements(eleid, :);
+      data(nanvflag(i)) = mean(data(nodeid(:)) , 'omitnan');
+   end
+   nanvflag = find(isnan(data));
+   NNanv = length(nanvflag);
+	count = count + 1;
+end
+newdata = data;
+	
Index: issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/integrateOverDomain.m
===================================================================
--- issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/integrateOverDomain.m	(revision 26975)
+++ issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/integrateOverDomain.m	(revision 26975)
@@ -0,0 +1,32 @@
+function [intData, meanData, areas] = integrateOverDomain(md, data, masked, weights)
+% integrateOverDomain - integrating data over the whole domain
+%
+%   intData: integral of the data over each element
+%   meanData: intData/areas
+%   areas: areas of the domain
+if nargin < 4
+	weights = ones(size(data));
+	if nargin<3
+		masked = [];
+	end
+end
+
+% Set the area with masked=1 to nan
+data(masked) = nan;
+weights(masked) =nan;
+
+% get the mesh
+elements=md.mesh.elements;
+x=md.mesh.x;
+y=md.mesh.y;
+
+%compute areas;
+eleAreas=GetAreas(elements,x,y);
+
+% integrate nodal data to element
+eleData = 1/3*eleAreas.*(data(elements(:,1),:).*weights(elements(:,1),:) + data(elements(:,2),:).*weights(elements(:,2),:) + data(elements(:,3),:).*weights(elements(:,3),:));
+eleAreas = 1/3*eleAreas.*(weights(elements(:,1),:)+weights(elements(:,2),:)+weights(elements(:,3),:));
+
+intData = sum(eleData(:),'omitnan');
+areas = sum(eleAreas(:),'omitnan');
+meanData = intData / areas;
Index: issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/integrateOverTime.m
===================================================================
--- issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/integrateOverTime.m	(revision 26975)
+++ issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/integrateOverTime.m	(revision 26975)
@@ -0,0 +1,26 @@
+function [newData, intData] = integrateOverTime(tdata, data, time, Nintdt)
+% integrateOverDomain - integrating data over the whole domain
+%
+%	 Input:
+%		tdata:		time steps of data
+%		data:			time series of the data 
+%		time:			time steps to be used for integration
+%		Nintdt:		number of time steps in 'time' to be used for integration
+%
+%	 Return:
+%		newData:		data projected on to 'time'
+%		intData:		integrated 'newData'	
+
+dt = abs(time(2) -time(1));
+% patch the first N time steps
+intTime = [[-Nintdt:-1]*dt+time(1), time];
+% integration scheme
+integ = ones(1, Nintdt) * dt;
+
+% interpolate data from the time scheme 'tdata' to 'time'
+newData = interp1(double(tdata), double(data), intTime, 'linear', 'extrap'); 
+
+% integrated, only works for equidistance time steps for now
+intData = conv(integ, newData);
+intData = intData(Nintdt:end-Nintdt);
+newData = newData(Nintdt+1: end);
Index: issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/interpZeroPos.m
===================================================================
--- issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/interpZeroPos.m	(revision 26975)
+++ issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/interpZeroPos.m	(revision 26975)
@@ -0,0 +1,52 @@
+function X0 = interpZeroPos(X, Y)
+%interpZeroPos - find the x values for y=0 with line interpolation of the
+%                given X and Y coordinates.  
+%                   X: is a matrix of size mxn, each column represent a dimension
+%                   Y: is a matrix of size mxk, k is the number of date sets, and each 
+%                       data set must contains positive and negative values.
+%                 The return value:
+%                   x0: is a matrix of size kxn, the rows correspond to each column
+%                   of X and the columns are for each data set.
+%                
+%               The method is to find the two consective values y1 and y2 where
+%               Y changes the sign, and the corresponding X are x1 and x2.
+%                      x2*y1-x1*y2
+%                x0 = -------------
+%                        y1-y2
+% Author: Cheng Gong
+% Last modified: 2020-10-20
+[m, n] = size(X);
+[my, k] = size(Y);
+if m~=my
+    error('rows of X and Y are not the same!');
+end
+
+% find the place where Y changes signs, or y=0 at some points
+mask = (Y(1:m-1,:).*Y(2:m,:) <=0);
+[row, col] = find(mask);
+ind = row;
+
+% use only the first index from each column
+if length(col)>k
+    disp('Multiple candidates of the y=0 values are found in one data set! Be defalut, take the first y=0 from upstream.');
+    [col, ia, ~] = unique(col);
+    ind = row(ia);
+end
+
+if length(col)<k
+	disp('y does not contain both positive and negative values, find the closest value to 0');
+	[~,ind] = min(abs(Y(1:m-1,:)));
+	ind = ind';
+	col= [1:k]';
+	if (length(ind) ~= length(col))
+		error('Size of the row and col indicators are not consistent!');
+	end
+end
+
+a1 = zeros(1,k);
+for i =1:k
+    a1(i) =  Y(ind(i),col(i))./(Y(ind(i),col(i))-Y(ind(i)+1,col(i)));
+end
+
+a2 = 1- a1; 
+X0 = X(ind+1,:).*a1(:) + X(ind,:).*a2(:);
Index: issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/inverseHistograms.m
===================================================================
--- issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/inverseHistograms.m	(revision 26975)
+++ issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/inverseHistograms.m	(revision 26975)
@@ -0,0 +1,34 @@
+function [cdfa, q_Sinv] = inverseHistograms(mu, S, Sinv)
+%     Given a distribution mu compute its inverse quantile function
+%     Parameters
+%     ----------
+%     mu     : histogram
+%     S      : support of the histogram
+%     Sinv   : support of the quantile function
+%     Returns
+%     -------
+%     cdfa   : the cumulative distribution function and
+%     q_Sinv : the inverse quantile function of the distribution mu
+
+epsilon = 1e-14;
+A = (mu>epsilon);
+Sa = S(A);
+
+% cummulative sum
+cdf = cumtrapz(S, mu);
+cdfa = cdf(A);
+
+% set the first value to 0 and last value to 1
+cdfa(1) = 0;
+Sa(1) = 0;
+
+if cdfa(end) < 1
+    cdfa = [cdfa(:);1];
+    Sa = [Sa(:);S(end)];
+end
+
+[~, ind] = unique(cdfa);
+% linear interpolation
+q_Sinv = interp1(cdfa(ind), Sa(ind), Sinv);
+
+end
Index: issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/movingAverage.m
===================================================================
--- issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/movingAverage.m	(revision 26975)
+++ issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/movingAverage.m	(revision 26975)
@@ -0,0 +1,42 @@
+function outdata =  movingAverage(rawdata, varargin)
+%
+% rawdata is an MxN matrix, where the first M-1 rows are the raw data, the row M is the time
+
+options    = pairoptions(varargin{:});
+timeWindow	= getfieldvalue(options, 'time window', 30); % unit in days
+resamp = getfieldvalue(options, 'resample', 1); % resample the data set with the size of window
+
+% get data info {{{
+[M, N] = size(rawdata);
+if M < 2
+	error('The data has to be at least two rows!');
+end
+if N < 2
+	error('The data has to be at least two columns!');
+end
+
+% sort the time
+time = rawdata(M, :);
+[time, tI] = sort(time);
+data = rawdata(1:M-1, :);
+data = data(:, tI);
+
+dt = time(2) - time(1);
+if dt < 0
+	error('The time need to be in ascending order.');
+end
+windowSize = ceil(timeWindow/365/dt);
+%}}}
+%% filter {{{
+if resamp
+	ind = [1:timeWindow:N];
+else
+	ind = [1:N];
+end
+b = (1/windowSize)*ones(1,windowSize);
+a = 1;
+disp(['Applying filter of width ', num2str(timeWindow)]);
+smoothdata = filter(b,a,data,[],2);
+time = time(ind);
+outdata = [smoothdata(:,ind);time];
+%}}}
Index: issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/npsd.m
===================================================================
--- issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/npsd.m	(revision 26975)
+++ issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/npsd.m	(revision 26975)
@@ -0,0 +1,32 @@
+function [npsdx, psdx, freq] = npsd(x)
+%npsd - Normalized Power Spectral Density (using fft)
+%
+%   x: time series
+%
+%   psdx: Power Spectral Density(right half plane, value doubled)
+%   freq: frequency(non-negative half)
+%
+%   Author: Cheng Gong
+%   Date: 2021-08-17
+
+% length of data
+Nx = length(x);
+% fft
+xdft = fft(x);
+xdft = xdft(1:floor(Nx/2+1));
+% power specturm
+psdx = (1/(2*pi*Nx)) * abs(xdft).^2;
+%  In order to conserve the total power, multiply all frequencies that
+%  occur in both sets - the positive and negative frequencies — by a factor of 2
+psdx(2:end-1) = 2*psdx(2:end-1);
+freq = 0:(2*pi)/Nx:pi;
+% normalize(trapzoidal rule)
+intS = trapz(freq, psdx);
+if intS == 0
+    psdx = 1 +psdx;
+    intS = trapz(freq, psdx);
+end
+npsdx = psdx ./ intS;
+
+
+end
Index: issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/projectToFlowlines.m
===================================================================
--- issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/projectToFlowlines.m	(revision 26975)
+++ issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/projectToFlowlines.m	(revision 26975)
@@ -0,0 +1,16 @@
+% Function to project values on the mesh to a flowline
+%
+%	md			-	ISSM model with mesh
+%	pValue	-	data on the mesh
+%	fx			-	x coordinates of the flowline
+%	fy			-	y coordinates of the flowline
+%
+% Author: Cheng Gong
+% Last modified: 2021-01-27
+
+function valueC = projectToFlowlines(md, pValue, fx, fy)
+    temp = InterpFromMeshToMesh2d(md.mesh.elements, md.mesh.x, md.mesh.y,...
+        pValue, fx, fy);
+    temp = diag(temp); % not very efficient, but works
+    valueC = temp(:)';
+end
Index: issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/psd.m
===================================================================
--- issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/psd.m	(revision 26975)
+++ issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/psd.m	(revision 26975)
@@ -0,0 +1,24 @@
+function [psdx, freq] = psd(x)
+%psd - Power Spectral Density (using fft)
+%
+%   x: time series
+%   
+%   psdx: Power Spectral Density(right half plane, value doubled)
+%   freq: frequency(non-negative half)
+%
+%   Author: Cheng Gong
+%   Date: 2021-08-17
+
+% length of data
+Nx = length(x);
+% fft
+xdft = fft(x);
+xdft = xdft(1:floor(Nx/2+1));
+% power specturm
+psdx = (1/(2*pi*Nx)) * abs(xdft).^2;
+%  In order to conserve the total power, multiply all frequencies that 
+%  occur in both sets - the positive and negative frequencies — by a factor of 2
+psdx(2:end-1) = 2*psdx(2:end-1);
+freq = 0:(2*pi)/Nx:pi;
+
+end
Index: issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/rescalegradientNan.m
===================================================================
--- issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/rescalegradientNan.m	(revision 26975)
+++ issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/rescalegradientNan.m	(revision 26975)
@@ -0,0 +1,7 @@
+% function to replace Nan by 0 in J, then rescale with the mass matrix
+%
+function J = rescalegradientNan(md, unscaledJ)
+
+nanflag = isnan(unscaledJ);
+unscaledJ(nanflag) = 0;
+J = rescalegradient(md, unscaledJ);
Index: issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/wfDistance.m
===================================================================
--- issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/wfDistance.m	(revision 26975)
+++ issm/trunk-jpl/src/m/contrib/chenggong/dataprocessing/wfDistance.m	(revision 26975)
@@ -0,0 +1,26 @@
+function w2 = wfDistance(x, y)
+%wfDistance - compute Wasserstein-Fourier Distance, according to Cazelles, E.,
+%         Robert, A. & Tobar, F. The Wasserstein-Fourier Distance for 
+%         Stationary Time Series. Ieee T Signal Proces 69, 709–721 (2019).
+%
+%   x: first time series
+%   y: second time series
+%
+%   w2:  Wasserstein-Fourier Distance, or W2 distance of sx and sy
+%
+%   Author: Cheng Gong
+%   Date: 2021-08-18
+N = length(x);
+Sinv=linspace(0,1,N);
+
+[csx, ~, freqx] = npsd(x);
+[csy, ~, freqy] = npsd(y);
+
+% inverse cummulative function
+[~, qx] = inverseHistograms(csx, freqx, Sinv);
+[~, qy] = inverseHistograms(csy, freqy, Sinv);
+
+%
+w2 = sqrt(trapz(Sinv ,(qx-qy).^2));
+
+end
Index: issm/trunk-jpl/src/m/contrib/chenggong/expxy2shpll.m
===================================================================
--- issm/trunk-jpl/src/m/contrib/chenggong/expxy2shpll.m	(revision 26975)
+++ issm/trunk-jpl/src/m/contrib/chenggong/expxy2shpll.m	(revision 26975)
@@ -0,0 +1,69 @@
+function expxy2shpll(expfilename,shpfilename,geometry)
+%expxy2shpll
+%  Convert .exp to .shp file
+%
+%   Usage:
+%			expxy2shpll('glacier.exp', 'glacier.shp')
+%			expxy2shpll('glacier.exp', 'glacier.shp', geometry)
+%
+%		geometry (optional)-	'MultiPoint' : point clouds
+%									'Point' : single point
+%									'Line' : two points line
+%									'polygon' : multiple points
+%
+
+%check file extensions
+[pathstr,name,ext] = fileparts(shpfilename);
+if ~strcmp(ext,'.shp'),
+   error(['Shapefile ' shpfilename ' does not have an extension .shp']);
+end
+
+[pathstr,name,ext] = fileparts(expfilename);
+if ~strcmp(ext,'.exp'),
+   error(['Exp file ' expfilename ' does not have an extension .exp']);
+end
+
+shp=expread(expfilename);
+[lat, lon] = xy2ll(shp.x, shp.y, 1, 45, 70);
+shp.x = lon;
+shp.y = lat;
+
+%initialize number of profile
+count=1;
+
+contours=struct([]);
+for i=1:length(shp),
+   if nargin < 3
+
+      %TEMP
+      %if contains(shp(i).name,'_pointcloud');
+      %  continue;
+      %end
+
+      if length(shp(i).x) == 0
+         continue;
+      elseif contains(shp(i).name,'_pointcloud');
+         geometry = 'MultiPoint';
+         shp(i).name = erase(shp(i).name,'_pointcloud');
+      elseif length(shp(i).x) == 1
+         geometry = 'Point';
+      elseif length(shp(i).x) < 3
+         geometry = 'Line';
+      else
+         if (shp(i).x(end)==shp(i).x(1) && shp(i).y(end)==shp(i).y(1)),
+            geometry = 'Polygon';
+         else
+            geometry = 'Line';
+         end
+      end
+   end
+   contours(count).Geometry=geometry;
+   contours(count).id=i;
+   contours(count).Name=shp(i).name;
+   contours(count).X=shp(i).x;
+   contours(count).Y=shp(i).y;
+   count = count+1;
+end
+
+%Make sure it is one single geometry otherwise it will yell at you
+shapewrite(contours,shpfilename);
Index: issm/trunk-jpl/src/m/contrib/chenggong/interpFromMEaSUREsGeotiff.m
===================================================================
--- issm/trunk-jpl/src/m/contrib/chenggong/interpFromMEaSUREsGeotiff.m	(revision 26973)
+++ issm/trunk-jpl/src/m/contrib/chenggong/interpFromMEaSUREsGeotiff.m	(revision 26975)
@@ -1,9 +1,33 @@
-% This function calls src/m/contrib/morlighem/modeldata/interpFromGeotiff.m for multiple times to load all avaliable 
-% tif data in  /totten_1/ModelData/Greenland/VelMEaSUREs/Jakobshavn_2008_2021/ within the given time period (in decimal years)
-% For some reason, each .tif file in this folder contains two sets of data, only the first dataset is useful
+function dataout = interpFromMEaSUREsGeotiff(X,Y,Tstart,Tend,varargin)
+%interpFromMEaSUREsGeotiff: 
+%	This function calls src/m/contrib/morlighem/modeldata/interpFromGeotiff.m for multiple times to load all avaliable 
+%	tif data in  /totten_1/ModelData/Greenland/VelMEaSUREs/Jakobshavn_2008_2021/ within the given time period (in decimal years)
+%	For some reason, each .tif file in this folder contains two sets of data, only the first dataset is useful
+%
+%   Usage:
+%		 dataout = interpFromMEaSUREsGeotiff(X,Y,Tstart,Tend, varargin)
+%
+%	X, Y are the coordinates of the mesh 
+%	Tstart and Tend decimal year of the start and end time
+%
+%   Example:
+%			obsData = interpFromMEaSUREsGeotiff(md.mesh.x,md.mesh.y, tstart, tend);
+%
+%   Options:
+%      - 'glacier':  which glacier to look for
+options    = pairoptions(varargin{:});
+glacier    = getfieldvalue(options,'glacier','Jakobshavn');
 
-function dataout = interpFromMEaSUREsGeotiff(X,Y,Tstart,Tend)
-
-foldername = '/totten_1/ModelData/Greenland/VelMEaSUREs/Jakobshavn_2008_2021/';
+if strcmp(glacier, 'Jakobshavn')
+	foldername = '/totten_1/ModelData/Greenland/VelMEaSUREs/Jakobshavn_2008_2021/';
+elseif strcmp(glacier, 'Kangerlussuaq')
+	foldername = '/totten_1/ModelData/Greenland/VelMEaSUREs/Kangerlussuaq_2006_2021/';
+elseif strcmp(glacier, 'Store')
+	foldername = '/totten_1/ModelData/Greenland/VelMEaSUREs/Store_2008_2021/';
+elseif strcmp(glacier, 'Rink')
+	foldername = '/totten_1/ModelData/Greenland/VelMEaSUREs/Rink_2008_2022/';
+else
+	error(['The velocity data for ', glacier, ' is not available, please download from NSIDC first.']);
+end
 
 % get the time info from file names
Index: issm/trunk-jpl/src/m/contrib/chenggong/modeldata/interpISMIP6Temp.m
===================================================================
--- issm/trunk-jpl/src/m/contrib/chenggong/modeldata/interpISMIP6Temp.m	(revision 26975)
+++ issm/trunk-jpl/src/m/contrib/chenggong/modeldata/interpISMIP6Temp.m	(revision 26975)
@@ -0,0 +1,15 @@
+function temp = interpISMIP6Temp(X, Y)
+%interpISMIP6Temp - interpolate land ice basal temperature to the given mesh
+%
+%	X and Y are the coordinates of the mesh
+%
+%  Author: Cheng Gong
+%  Last modified: 2021-12-06
+
+filename = '/totten_1/ModelData/Greenland/ISMIP6/GreenlandISMIP6-Morlighem-2020-10-01.nc';
+
+x = ncread(filename, 'x');
+y = ncread(filename, 'y');
+tb = ncread(filename, 'tb');
+
+temp = InterpFromGrid(double(x), double(y), double(tb'), X, Y);
Index: issm/trunk-jpl/src/m/contrib/chenggong/visualization/imageNonUni.m
===================================================================
--- issm/trunk-jpl/src/m/contrib/chenggong/visualization/imageNonUni.m	(revision 26975)
+++ issm/trunk-jpl/src/m/contrib/chenggong/visualization/imageNonUni.m	(revision 26975)
@@ -0,0 +1,23 @@
+function imageNonUni(varargin)
+%imageNonUni - draw an image with non uniform grid
+%
+% Author: Cheng Gong
+% Last modified: 2020-08-24
+x = varargin{1};
+y = varargin{2};
+C = varargin{3};
+
+nx = length(x);
+ny = length(y);
+[ncy, ncx] = size(C);
+
+if ((nx ~= ncx) || (ny~=ncy))
+    error('The size of x-y coordinates do not match the data!');
+end
+
+X = repmat(x(:)',ny , 1);
+Y = repmat(y(:), 1, nx);
+
+
+surf(X, Y, zeros(ny,nx), C, 'EdgeColor','None');
+view(2);
Index: issm/trunk-jpl/src/m/contrib/chenggong/visualization/plotCompareTransientFlowline.m
===================================================================
--- issm/trunk-jpl/src/m/contrib/chenggong/visualization/plotCompareTransientFlowline.m	(revision 26975)
+++ issm/trunk-jpl/src/m/contrib/chenggong/visualization/plotCompareTransientFlowline.m	(revision 26975)
@@ -0,0 +1,48 @@
+function plotCompareTransientFlowline(velList, titleList, md, flowline, TStart, TEnd, xl, yl, ca)
+%plotCompareTransientFlowline - to compare the solutions along a given
+%  flowine for the given time periods
+%
+%   velList:	a cell of transient velocities on the model mesh
+%   titleList: names of the velocity data
+%   md:			ISSM model
+%   flowline:	a flowline data structure contains flowline.x and flowline.y. 
+%					If possible, it should also contains flowline.Xmain, which is 
+%					the distance along the flowline from upstream to the terminus position.
+%   TStart:		a list of start time points
+%   TEnd:		a list of end time points
+%   ca:			caixs value in plot
+%   xl:			xlim value
+%   yl:			ylim value
+%
+%   Author: Cheng Gong
+%   Date: 2021-12-06
+
+N = length(velList);
+vel_flowline = cell(N,1);
+if ~isfield(flowline, 'Xmain')
+    flowline.Xmain = cumsum([0;sqrt((flowline.x(2:end)- flowline.x(1:end-1)).^2 + (flowline.y(2:end)- flowline.y(1:end-1)).^2)]')/1000;
+end
+
+for i = 1:N
+    vel_flowline{i} = InterpFromMeshToMesh2d(md.mesh.elements,md.mesh.x,md.mesh.y,velList{i},flowline.x,flowline.y);
+end
+
+for p = 1:N
+    subplot(N,1,p)
+    hold on
+    for i = 1:length(TStart)
+        overlay = [vel_flowline{p}(:,i), vel_flowline{p}(:,i)]';
+        imAlpha = ones(size(overlay));
+        imAlpha(isnan(overlay)) = 0;
+        imageNonUni(flowline.Xmain, [TStart(i), TEnd(i)], overlay);
+    end
+    title(titleList{p}, 'Interpreter', 'latex');
+    xlim(xl)
+    xlabel('x(km)')
+    ylabel('year')
+    ylim(yl);
+    caxis(ca)
+    colormap('jet')
+    h=colorbar;
+    title(h,'m/a')
+end
