Function grdflexure_write_transfer_function

src/potential/grdflexure.c:1136–1197

}

GMT_LOCAL int grdflexure_write_transfer_function (struct GMT_CTRL *GMT, struct GRDFLEXURE_CTRL *Ctrl, struct GRDFLEXURE_RHEOLOGY *R, struct GMT_OPTION *options) {
	/* Write a table with six segments (one each for Te = 1, 2, 5, 10, 20, 50, and 100 km).
	 * Each segment has leading columns of wavelength and wavenumber corresponding to wavelengths 1:5000 km.
	 * The next 12 columns has the chosen transfer function evaluated for times 1k, 2k, 5k, 10k, 20k, 50k, 100k, 200k, 500k, 1M, 2M, and 5M years.
	 * Each segment is written to a separate file. Obviously, if no -F or -M are given then all columns are the same since elastic */
	uint64_t k;
	int t, s, n_times, n_te;
	char file[GMT_LEN64] = {""};
	static char *FLX_response[7] = {"Airy", "Elastic", "Viscoelastic", "Firmoviscous (1 layer)", "Firmoviscous (2 layer)", "Viscous (1 layer)", "Viscous (2 layer)"};
	uint64_t dim[4] = {1, 0, 0, 0};
	double *kr, K[3], te[7] = {1.0, 2.0, 5.0, 10.0, 20.0, 50.0, 100.0};
	double times[12] = {1.0, 2.0, 5.0, 10.0, 20.0, 50.0, 100.0, 200.0, 500.0, 1000.0, 2000.0, 5000.0};	/* Times in kiloyears */
	struct GMT_DATASET *D = NULL;
	struct GMT_DATASEGMENT *S = NULL;
	struct GMT_DATASEGMENT_HIDDEN *SH = NULL;
	struct GMT_ARRAY T;

	gmt_M_memset (&T, 1, struct GMT_ARRAY);	/* Wipe clean the structure */

	R->relative = true;	/* Relative times are implicitly given */
	n_te = (R->mode > FLX_FV2) ? 1 : 7;	/* For purely viscous we don't need to loop over plate thickness */
	gmt_parse_array (GMT, 'T', "1/5000/1", &T, GMT_ARRAY_RANGE | GMT_ARRAY_UNIQUE, 0);	/* In km */
	gmt_create_array (GMT, 'T', &T, NULL, NULL);
	dim[GMT_ROW] = T.n;
	n_times = (Ctrl->F.active || Ctrl->M.active) ? 12 : 1;	/* No point repeating 12 identical results for the elastic case */
	dim[GMT_SEG] = n_te;
	dim[GMT_COL] = 2 + n_times;
	GMT_Report (GMT->parent, GMT_MSG_DEBUG, "Write transfer functions\n");
	kr = gmt_M_memory (GMT, NULL, T.n, double);
	for (k = 0; k < T.n; k++) kr[k] = 2.0 * M_PI / (T.array[k] * 1000.0);	/* Radial wavenumber in 1/m */

	if ((D = GMT_Create_Data (GMT->parent, GMT_IS_DATASET, GMT_IS_LINE, 0, dim, NULL, NULL, 0, 0, NULL)) == NULL)
		return GMT_RUNTIME_ERROR;
	GMT_Set_Comment (GMT->parent, GMT_IS_DATASET, GMT_COMMENT_IS_OPTION | GMT_COMMENT_IS_COMMAND, options, D);
	gmt_set_tableheader (GMT, GMT_OUT, true);


	for (s = 0; s < n_te; s++) {
		S = D->table[0]->segment[s];
		Ctrl->E.te[TE_INIT] = (n_te > 1) ? te[s] * 1000 : 0.0;	/* Te in meters, zero for viscous only */
		sprintf (file, "grdflexure_transfer_function_te_%3.3d_km.txt", irint (Ctrl->E.te[TE_INIT] * 0.001));
		SH = gmt_get_DS_hidden (S);
		SH->file[GMT_OUT] = strdup (file);
		gmt_M_memcpy (S->data[0], T.array, T.n, double);
		gmt_M_memcpy (S->data[1], kr, T.n, double);
		GMT_Report (GMT->parent, GMT_MSG_INFORMATION, "%s transfer function for Te = %g km written to %s\n", FLX_response[R->mode], Ctrl->E.te[TE_INIT] * 0.001, SH->file[GMT_OUT]);
		for (t = 0; t < n_times; t++) {	/* For each time step (i.e., at least once) */
			R->eval_time_yr = times[t] * 1000.0;	/* In years */
			R->setup (GMT, Ctrl, R);		/* Set up parameters */
			R->scale = 1.0;	/* We want these to go 0-1 only */
			for (k = 0; k < T.n; k++) {	/* Evaluate transfer functions */
				K[GMT_FFT_K_IS_KR] = kr[k];
				S->data[t+2][k] = R->transfer (K, R);
			}
		}
	}
	gmt_free_array (GMT, &T);
	gmt_M_free (GMT, kr);