onion.cpp @ 8f5b34a

ESS_GUIESS_GUI_DocsESS_GUI_batch_fittingESS_GUI_bumps_abstractionESS_GUI_iss1116ESS_GUI_iss879ESS_GUI_iss959ESS_GUI_openclESS_GUI_orderingESS_GUI_sync_sascalccostrafo411magnetic_scattrelease-4.1.1release-4.1.2release-4.2.2release_4.0.1ticket-1009ticket-1094-headlessticket-1242-2d-resolutionticket-1243ticket-1249ticket885unittest-saveload

Last change on this file since 8f5b34a was 34c2649, checked in by Mathieu Doucet <doucetm@…>, 13 years ago
Re #4 Fixed warnings
Property mode set to `100644`
File size: 14.1 KB

Rev	Line
[e096270]	1	/**
	2	This software was developed by the University of Tennessee as part of the
	3	Distributed Data Analysis of Neutron Scattering Experiments (DANSE)
	4	project funded by the US National Science Foundation.
	5
	6	If you use DANSE applications to do scientific research that leads to
	7	publication, we ask that you acknowledge the use of the software with the
	8	following sentence:
	9
	10	"This work benefited from DANSE software developed under NSF award DMR-0520547."
	11
	12	copyright 2008, University of Tennessee
	13	*/
	14
	15	/**
	16	* Scattering model classes
	17	* The classes use the IGOR library found in
	18	* sansmodels/src/libigor
	19	*
	20	*/
	21
	22	#include <math.h>
	23	#include "models.hh"
	24	#include "parameters.hh"
	25	#include <stdio.h>
	26	using namespace std;
	27
	28	extern "C" {
	29	#include "onion.h"
	30	}
	31
	32	OnionModel :: OnionModel() {
	33	n_shells = Parameter(1.0);
	34	scale = Parameter(1.0);
[339ce67]	35	rad_core0 = Parameter(200.0);
	36	sld_core0 = Parameter(1e-06);
[e096270]	37	sld_solv = Parameter(6.4e-06);
	38	background = Parameter(0.0);
	39
	40
	41	sld_out_shell1 = Parameter(1.0e-06);
	42	sld_out_shell2 = Parameter(1.0e-06);
	43	sld_out_shell3 = Parameter(1.0e-06);
	44	sld_out_shell4 = Parameter(1.0e-06);
	45	sld_out_shell5 = Parameter(1.0e-06);
	46	sld_out_shell6 = Parameter(1.0e-06);
	47	sld_out_shell7 = Parameter(1.0e-06);
	48	sld_out_shell8 = Parameter(1.0e-06);
	49	sld_out_shell9 = Parameter(1.0e-06);
	50	sld_out_shell10 = Parameter(1.0e-06);
	51
	52
	53	sld_in_shell1 = Parameter(2.3e-06);
	54	sld_in_shell2 = Parameter(2.6e-06);
	55	sld_in_shell3 = Parameter(2.9e-06);
	56	sld_in_shell4 = Parameter(3.2e-06);
	57	sld_in_shell5 = Parameter(3.5e-06);
	58	sld_in_shell6 = Parameter(3.8e-06);
	59	sld_in_shell7 = Parameter(4.1e-06);
	60	sld_in_shell8 = Parameter(4.4e-06);
	61	sld_in_shell9 = Parameter(4.7e-06);
	62	sld_in_shell10 = Parameter(5.0e-06);
	63
	64
	65	A_shell1 = Parameter(1.0);
	66	A_shell2 = Parameter(1.0);
	67	A_shell3 = Parameter(1.0);
	68	A_shell4 = Parameter(1.0);
	69	A_shell5 = Parameter(1.0);
	70	A_shell6 = Parameter(1.0);
	71	A_shell7 = Parameter(1.0);
	72	A_shell8 = Parameter(1.0);
	73	A_shell9 = Parameter(1.0);
	74	A_shell10 = Parameter(1.0);
	75
	76
	77	thick_shell1 = Parameter(50.0);
	78	thick_shell2 = Parameter(50.0);
	79	thick_shell3 = Parameter(50.0);
	80	thick_shell4 = Parameter(50.0);
	81	thick_shell5 = Parameter(50.0);
	82	thick_shell6 = Parameter(50.0);
	83	thick_shell7 = Parameter(50.0);
	84	thick_shell8 = Parameter(50.0);
	85	thick_shell9 = Parameter(50.0);
	86	thick_shell10 = Parameter(50.0);
	87
	88
	89	func_shell1 = Parameter(2);
	90	func_shell2 = Parameter(2);
	91	func_shell3 = Parameter(2);
	92	func_shell4 = Parameter(2);
	93	func_shell5 = Parameter(2);
	94	func_shell6 = Parameter(2);
	95	func_shell7 = Parameter(2);
	96	func_shell8 = Parameter(2);
	97	func_shell9 = Parameter(2);
	98	func_shell10 = Parameter(2);
	99
	100	}
	101
	102	/**
	103	* Function to evaluate 1D scattering function
	104	* The NIST IGOR library is used for the actual calculation.
	105	* @param q: q-value
	106	* @return: function value
	107	*/
	108	double OnionModel :: operator()(double q) {
	109	double dp[56];
	110	// Fill parameter array for IGOR library
	111	// Add the background after averaging
	112	dp[0] = n_shells();
	113	dp[1] = scale();
[339ce67]	114	dp[2] = rad_core0();
	115	dp[3] = sld_core0();
[e096270]	116	dp[4] = sld_solv();
	117	dp[5] = 0.0;
	118
	119	dp[6] = sld_out_shell1();
	120	dp[7] = sld_out_shell2();
	121	dp[8] = sld_out_shell3();
	122	dp[9] = sld_out_shell4();
	123	dp[10] = sld_out_shell5();
	124	dp[11] = sld_out_shell6();
	125	dp[12] = sld_out_shell7();
	126	dp[13] = sld_out_shell8();
	127	dp[14] = sld_out_shell9();
	128	dp[15] = sld_out_shell10();
	129
	130	dp[16] = sld_in_shell1();
	131	dp[17] = sld_in_shell2();
	132	dp[18] = sld_in_shell3();
	133	dp[19] = sld_in_shell4();
	134	dp[20] = sld_in_shell5();
	135	dp[21] = sld_in_shell6();
	136	dp[22] = sld_in_shell7();
	137	dp[23] = sld_in_shell8();
	138	dp[24] = sld_in_shell9();
	139	dp[25] = sld_in_shell10();
	140
	141	dp[26] = A_shell1();
	142	dp[27] = A_shell2();
	143	dp[28] = A_shell3();
	144	dp[29] = A_shell4();
	145	dp[30] = A_shell5();
	146	dp[31] = A_shell6();
	147	dp[32] = A_shell7();
	148	dp[33] = A_shell8();
	149	dp[34] = A_shell9();
	150	dp[35] = A_shell10();
	151
	152	dp[36] = thick_shell1();
	153	dp[37] = thick_shell2();
	154	dp[38] = thick_shell3();
	155	dp[39] = thick_shell4();
	156	dp[40] = thick_shell5();
	157	dp[41] = thick_shell6();
	158	dp[42] = thick_shell7();
	159	dp[43] = thick_shell8();
	160	dp[44] = thick_shell9();
	161	dp[45] = thick_shell10();
	162
	163	dp[46] = func_shell1();
	164	dp[47] = func_shell2();
	165	dp[48] = func_shell3();
	166	dp[49] = func_shell4();
	167	dp[50] = func_shell5();
	168	dp[51] = func_shell6();
	169	dp[52] = func_shell7();
	170	dp[53] = func_shell8();
	171	dp[54] = func_shell9();
	172	dp[55] = func_shell10();
	173
	174
	175	// Get the dispersion points for the radius
	176	vector<WeightPoint> weights_rad;
[339ce67]	177	rad_core0.get_weights(weights_rad);
[e096270]	178
	179	// Get the dispersion points for the thick 1
	180	vector<WeightPoint> weights_s1;
	181	thick_shell1.get_weights(weights_s1);
	182
	183	// Get the dispersion points for the thick 2
	184	vector<WeightPoint> weights_s2;
	185	thick_shell2.get_weights(weights_s2);
	186
	187	// Get the dispersion points for the thick 3
	188	vector<WeightPoint> weights_s3;
	189	thick_shell3.get_weights(weights_s3);
	190
	191	// Get the dispersion points for the thick 4
	192	vector<WeightPoint> weights_s4;
	193	thick_shell4.get_weights(weights_s4);
	194
	195	// Get the dispersion points for the thick 5
	196	vector<WeightPoint> weights_s5;
	197	thick_shell5.get_weights(weights_s5);
	198
	199	// Get the dispersion points for the thick 6
	200	vector<WeightPoint> weights_s6;
	201	thick_shell6.get_weights(weights_s6);
	202
	203	// Get the dispersion points for the thick 7
	204	vector<WeightPoint> weights_s7;
	205	thick_shell7.get_weights(weights_s7);
	206
	207	// Get the dispersion points for the thick 8
	208	vector<WeightPoint> weights_s8;
	209	thick_shell8.get_weights(weights_s8);
	210	// Get the dispersion points for the thick 9
	211	vector<WeightPoint> weights_s9;
	212	thick_shell9.get_weights(weights_s9);
	213
	214	// Get the dispersion points for the thick 10
	215	vector<WeightPoint> weights_s10;
	216	thick_shell10.get_weights(weights_s10);
	217
	218
	219	// Perform the computation, with all weight points
	220	double sum = 0.0;
	221	double norm = 0.0;
	222	double vol = 0.0;
	223
	224	// Loop over radius weight points
[34c2649]	225	for(size_t i=0; i<weights_rad.size(); i++) {
[e096270]	226	dp[2] = weights_rad[i].value;
	227	// Loop over radius weight points
[34c2649]	228	for(size_t j=0; j<weights_s1.size(); j++) {
[e096270]	229	dp[36] = weights_s1[j].value;
	230	// Loop over radius weight points
[34c2649]	231	for(size_t k=0; k<weights_s2.size(); k++) {
[e096270]	232	dp[37] = weights_s2[k].value;
	233	// Loop over radius weight points
[34c2649]	234	for(size_t l=0; l<weights_s3.size(); l++) {
[e096270]	235	dp[38] = weights_s3[l].value;
	236	// Loop over radius weight points
[34c2649]	237	for(size_t m=0; m<weights_s4.size(); m++) {
[e096270]	238	dp[39] = weights_s4[m].value;
[34c2649]	239	for(size_t n=0; n<weights_s5.size(); n++) {
[e096270]	240	dp[40] = weights_s5[n].value;
[34c2649]	241	for(size_t o=0; o<weights_s6.size(); o++) {
[e096270]	242	dp[41] = weights_s6[o].value;
[34c2649]	243	for(size_t p=0; p<weights_s7.size(); p++) {
[e096270]	244	dp[42] = weights_s7[p].value;
[34c2649]	245	for(size_t t=0; t<weights_s8.size(); t++) {
[e096270]	246	dp[43] = weights_s8[t].value;
[34c2649]	247	for(size_t r=0; r<weights_s9.size(); r++) {
[e096270]	248	dp[44] = weights_s9[r].value;
[34c2649]	249	for(size_t s=0; s<weights_s10.size(); s++) {
[e096270]	250	dp[45] = weights_s10[s].value;
	251	//Un-normalize Shells by volume
	252	sum += weights_rad[i].weightweights_s1[j].weightweights_s2[k].weightweights_s3[l].weightweights_s4[m].weight
	253	weights_s5[n].weightweights_s6[o].weightweights_s7[p].weightweights_s8[t].weight
	254	weights_s9[r].weightweights_s10[s].weight
	255	* so_kernel(dp,q) * pow((weights_rad[i].value+weights_s1[j].value+weights_s2[k].value+weights_s3[l].value+weights_s4[m].value
	256	+weights_s5[n].value+weights_s6[o].value+weights_s7[p].value+weights_s8[t].value
	257	+weights_s9[r].value+weights_s10[s].value),3.0);
	258	//Find average volume
	259	vol += weights_rad[i].weightweights_s1[j].weightweights_s2[k].weightweights_s3[l].weightweights_s4[m].weight
	260	weights_s5[n].weightweights_s6[o].weightweights_s7[p].weightweights_s8[t].weight
	261	weights_s9[r].weightweights_s10[s].weight
	262	* pow((weights_rad[i].value+weights_s1[j].value+weights_s2[k].value+weights_s3[l].value+weights_s4[m].value
	263	+weights_s5[n].value+weights_s6[o].value+weights_s7[p].value+weights_s8[t].value
	264	+weights_s9[r].value+weights_s10[s].value),3.0);
	265	norm += weights_rad[i].weightweights_s1[j].weightweights_s2[k].weightweights_s3[l].weightweights_s4[m].weight
	266	weights_s5[n].weightweights_s6[o].weightweights_s7[p].weightweights_s8[t].weight
	267	weights_s9[r].weightweights_s10[s].weight;
	268	}
	269	}
	270	}
	271	}
	272	}
	273	}
	274	}
	275	}
	276	}
	277	}
	278	}
	279
	280	if (vol != 0.0 && norm != 0.0) {
	281	//Re-normalize by avg volume
	282	sum = sum/(vol/norm);}
	283
	284	return sum/norm + background();
	285	}
	286
	287	/**
	288	* Function to evaluate 2D scattering function
	289	* @param q_x: value of Q along x
	290	* @param q_y: value of Q along y
	291	* @return: function value
	292	*/
	293	double OnionModel :: operator()(double qx, double qy) {
	294	double q = sqrt(qxqx + qyqy);
	295	return (*this).operator()(q);
	296	}
	297
	298	/**
	299	* Function to evaluate 2D scattering function
	300	* @param pars: parameters of the sphere
	301	* @param q: q-value
	302	* @param phi: angle phi
	303	* @return: function value
	304	*/
	305	double OnionModel :: evaluate_rphi(double q, double phi) {
	306	return (*this).operator()(q);
	307	}
	308
	309	/**
	310	* Function to calculate effective radius
	311	* @return: effective radius value
	312	*/
	313	double OnionModel :: calculate_ER() {
	314	OnionParameters dp;
[339ce67]	315	dp.rad_core0 = rad_core0();
[e096270]	316	dp.thick_shell1 = thick_shell1();
	317	dp.thick_shell2 = thick_shell2();
	318	dp.thick_shell3 = thick_shell3();
	319	dp.thick_shell4 = thick_shell4();
	320	dp.thick_shell5 = thick_shell5();
	321	dp.thick_shell6 = thick_shell6();
	322	dp.thick_shell7 = thick_shell7();
	323	dp.thick_shell8 = thick_shell8();
	324	dp.thick_shell9 = thick_shell9();
	325	dp.thick_shell10 = thick_shell10();
	326
	327
	328	double rad_out = 0.0;
	329	// Perform the computation, with all weight points
	330	double sum = 0.0;
	331	double norm = 0.0;
	332
	333	// Get the dispersion points for the radius
	334	vector<WeightPoint> weights_rad;
[339ce67]	335	rad_core0.get_weights(weights_rad);
[e096270]	336
	337	// Get the dispersion points for the thick 1
	338	vector<WeightPoint> weights_s1;
	339	thick_shell1.get_weights(weights_s1);
	340
	341	// Get the dispersion points for the thick 2
	342	vector<WeightPoint> weights_s2;
	343	thick_shell2.get_weights(weights_s2);
	344
	345	// Get the dispersion points for the thick 3
	346	vector<WeightPoint> weights_s3;
	347	thick_shell3.get_weights(weights_s3);
	348
	349	// Get the dispersion points for the thick 4
	350	vector<WeightPoint> weights_s4;
	351	thick_shell4.get_weights(weights_s4);
	352	// Get the dispersion points for the thick 5
	353	vector<WeightPoint> weights_s5;
	354	thick_shell5.get_weights(weights_s5);
	355
	356	// Get the dispersion points for the thick 6
	357	vector<WeightPoint> weights_s6;
	358	thick_shell6.get_weights(weights_s6);
	359
	360	// Get the dispersion points for the thick 7
	361	vector<WeightPoint> weights_s7;
	362	thick_shell7.get_weights(weights_s7);
	363
	364	// Get the dispersion points for the thick 8
	365	vector<WeightPoint> weights_s8;
	366	thick_shell8.get_weights(weights_s8);
	367	// Get the dispersion points for the thick 9
	368	vector<WeightPoint> weights_s9;
	369	thick_shell9.get_weights(weights_s9);
	370
	371	// Get the dispersion points for the thick 10
	372	vector<WeightPoint> weights_s10;
	373	thick_shell10.get_weights(weights_s10);
	374
	375
	376	// Loop over radius weight points
[34c2649]	377	for(size_t i=0; i<weights_rad.size(); i++) {
[339ce67]	378	dp.rad_core0 = weights_rad[i].value;
[e096270]	379	// Loop over radius weight points
[34c2649]	380	for(size_t j=0; j<weights_s1.size(); j++) {
[e096270]	381	dp.thick_shell1 = weights_s1[j].value;
	382	// Loop over radius weight points
[34c2649]	383	for(size_t k=0; k<weights_s2.size(); k++) {
[e096270]	384	dp.thick_shell2 = weights_s2[k].value;
	385	// Loop over radius weight points
[34c2649]	386	for(size_t l=0; l<weights_s3.size(); l++) {
[e096270]	387	dp.thick_shell3 = weights_s3[l].value;
	388	// Loop over radius weight points
[34c2649]	389	for(size_t m=0; m<weights_s4.size(); m++) {
[e096270]	390	dp.thick_shell4 = weights_s4[m].value;
	391	// Loop over radius weight points
[34c2649]	392	for(size_t n=0; j<weights_s5.size(); n++) {
[e096270]	393	dp.thick_shell5 = weights_s5[n].value;
	394	// Loop over radius weight points
[34c2649]	395	for(size_t o=0; k<weights_s6.size(); o++) {
[e096270]	396	dp.thick_shell6 = weights_s6[o].value;
	397	// Loop over radius weight points
[34c2649]	398	for(size_t p=0; l<weights_s7.size(); p++) {
[e096270]	399	dp.thick_shell7 = weights_s7[p].value;
	400	// Loop over radius weight points
[34c2649]	401	for(size_t t=0; m<weights_s8.size(); t++) {
[e096270]	402	dp.thick_shell8 = weights_s8[t].value;
	403	// Loop over radius weight points
[34c2649]	404	for(size_t r=0; l<weights_s9.size(); r++) {
[e096270]	405	dp.thick_shell8 = weights_s9[r].value;
	406	// Loop over radius weight points
[34c2649]	407	for(size_t s=0; m<weights_s10.size(); s++) {
[e096270]	408	dp.thick_shell10 = weights_s10[s].value;
	409	//Un-normalize FourShell by volume
	410	sum += weights_rad[i].weightweights_s1[j].weightweights_s2[k].weightweights_s3[l].weightweights_s4[m].weight
	411	weights_s5[n].weightweights_s6[o].weightweights_s7[p].weightweights_s8[t].weight
	412	weights_s9[r].weightweights_s10[s].weight
[339ce67]	413	* (dp.rad_core0+dp.thick_shell1+dp.thick_shell2+dp.thick_shell3+dp.thick_shell4+dp.thick_shell5
[e096270]	414	+dp.thick_shell6+dp.thick_shell7+dp.thick_shell8+dp.thick_shell9+dp.thick_shell10);
	415	norm += weights_rad[i].weightweights_s1[j].weightweights_s2[k].weight*weights_s3[l].weight
	416	weights_s4[m].weightweights_s5[n].weightweights_s6[o].weightweights_s7[p].weight
	417	weights_s8[t].weightweights_s9[r].weight*weights_s10[s].weight;
	418	}
	419	}
	420	}
	421	}
	422	}
	423	}
	424	}
	425	}
	426	}
	427	}
	428	}
	429
	430	if (norm != 0){
	431	//return the averaged value
	432	rad_out = sum/norm;}
	433	else{
	434	//return normal value
[339ce67]	435	rad_out = dp.rad_core0+dp.thick_shell1+dp.thick_shell2+dp.thick_shell3+dp.thick_shell4
[e096270]	436	+dp.thick_shell5+dp.thick_shell6+dp.thick_shell7+dp.thick_shell8+dp.thick_shell9+dp.thick_shell10;}
	437	return rad_out;
	438	}

Note: See TracBrowser for help on using the repository browser.

SasView

source: sasview/sansmodels/src/sans/models/c_models/onion.cpp @ 8f5b34a

Download in other formats: