Context Navigation

Back to Ticket #2156

Ticket #2156: COLETTE_GET_SPECIAL.F.TXT

File COLETTE_GET_SPECIAL.F.TXT, 24.5 KB (added by Nick Draper, 10 years ago)

Line
1	c===============================================================================
2	c COLETTE - LOQ data analysis development program - R.Osborn Nov. 1985
3	c
4	c Produce data in the Special form in workspace i_graph for subsequent
5	c plotting
6	c
7	integer*4 function get_special (n)
8
9	external cli$_absent
10	character*80 text
11	integer*4 option, type_of_plot, subtract_wkspc
12	logical*1 first_time /.TRUE./, hist
13
14	include 'COLETTE_SOURCES:COMSMG.INC'
15	include 'COLETTE_COMLOQ:'
16	COMMON/DSPECIAL/SP_E,SP_F,SP_G,SP_H,SP_CX,
17	> SP_A,SP_B,SP_C,SP_D,SP_CY
18	data SP_E,SP_F,SP_G,SP_H,SP_CX/0.,0.,1.,0.,1./,
19	> SP_A,SP_B,SP_C,SP_D,SP_CY/0.,0.,4.,1.,1./
20	c
21	c==== RKH 29/10/92
22	C==== pull in something to tell whether we are in POINT mode or not !
23	c==== will program this for all basic plotting options - should not
24	c==== be relevant for histogram mode !
25	c==== note COLETTE generates binned data with histogram x boundaries,
26	c==== but these do not behave as proper histograms as the Y values are strictly
27	c==== funtions NOT histograms !
28	c==== on writing out to file we take mid-points, so OLD'ed data really are
29	C==== point mode ! OLD by default takes mid-points again to generate more
30	C==== histogram x bins, this is OK for dq=constant, but not for log or
31	c==== irregular bins. Hence you can use OLD/POINT, then TOG MODE in GENIE
32	c==== after which D/S ought mow to work properly.
33	c
34	c==== include 'OGENIE_SOURCES:GKSPLOT.CMN'
35	c==== this has equivalent of:
36	common/drtype/type_of_plot, hist
37
38	c===============================================================================
39	c Output help screen for special plots
40
41	if (first_time) then
42	status = smg$create_virtual_display (12, 70, vdid8, smg$m_border)
43	status = smg$put_chars
44	> (vdid8, 'Option nos. of Special Plots of IvQ', 1, 4,, smg$m_underline)
45	status = smg$put_chars
46	> (vdid8, '0: EXIT THIS MENU ', 2, 4)
47	status = smg$put_chars
48	> (vdid8, '1: Guinier (spheres) - Ln(I) v Q**2', 3, 4)
49	status = smg$put_chars
50	> (vdid8, '2: Guinier (rods) - Ln(IQ) v Q**2', 4, 4)
51	status = smg$put_chars
52	> (vdid8, '3: Guinier (sheets) - Ln(IQ2) v Q2', 5, 4)
53	status = smg$put_chars
54	> (vdid8, '4: Zimm - 1/I v Q**2', 6, 4)
55	status = smg$put_chars
56	> (vdid8, '5: Kratky - IQ**2 v Q', 7, 4)
57	status = smg$put_chars
58	> (vdid8, '6: Debye-Bueche - 1/SqrtI v Q**2', 8, 4)
59	status = smg$put_chars
60	> (vdid8, '7: Log - Log Ln(I) v Ln(Q)', 9, 4)
61	status = smg$put_chars
62	> (vdid8, '8: Porod IQ**4 v Q',10, 4)
63	status = smg$put_chars
64	>(vdid8,'9: General (QA)(IB)xLn(QCxID) v similar',11, 4)
65	status = smg$put_chars
66	>(vdid8,'add 100 for true histogram; 91 re-uses last 9 ', 12, 4)
67	first_time = .FALSE.
68	end if
69
70	status = cli$get_value ('SPECIAL', text, len)
71
72	if (status .eq. %loc(cli$_absent)) then
73	status = smg$paste_virtual_display (vdid8, pbid, 7, 6)
74	call prompt ('Give option no ==>', text, len)
75	status = smg$unpaste_virtual_display (vdid8,pbid)
76	end if
77
78	if (len .eq. 0) then
79	get_special = status_not_ok
80	call replace_prompt
81	return
82	end if
83
84	read (text(1:len), '(i)', err=100) option
85	if (option.eq.0) then
86	get_special=status_not_ok
87	call replace_prompt
88	return
89	endif
90	c===============================================================================
91	c Subtract background level if required
92	* Modified to allow background to be given in a workspace, SMK, 02-11-92
93	* Modified to allow background to be given as command line qualifier, SMK
94
95	status = cli$get_value ('BACKGROUND', text, len)
96	if (status .eq. %loc(cli$_absent)) then
97	call prompt ('Background level ==> ', text, len)
98	end if
99
100	subtract_wkspc=0
101	* IF NOTHING GIVEN THEN ASSUME BACKGROUND IS ZERO
102	if (len.eq.0) then
103	bkgd=0.0
104	* OTHERWISE, CHECK TO SEE IF BACKGROUND IS CONTAINED IN A WORKSPACE
105	else if (text(1:1).eq.'W') then
106	if (len.eq.1) goto 100
107	read (text(2:len),fmt=*,err=100) j_graph
108	if ( (j_graph.le.0).or.
109	& ((j_graph.ge.15.).and.(j_graph.le.24)).or.
110	& ((j_graph.ge.61.).and.(j_graph.le.69)).or.
111	& (j_graph.gt.used_streams) ) then
112	call error('ERROR: That is a reserved workspace!')
113	goto 102
114	end if
115	if (npt(j_graph).le.0) then
116	call error('ERROR: That workspace is empty!')
117	goto 102
118	end if
119	offset_j_graph=len_streams*(j_graph-1)
120	subtract_wkspc=1
121	* IF BACKGROUND VALUE GIVEN
122	else
123	read (text(1:len),fmt=*,err=100) bkgd
124	end if
125
126	npt(i_graph) = npt(n)
127	offset_i_graph = len_streams * (i_graph-1)
128	offset_n = len_streams * (n - 1)
129	n1=offset_i_graph+1
130	n2=offset_i_graph+npt(i_graph)
131
132	c===============================================================================
133	c Perform workspace transformation according to option no.
134	c
135	c I(Q) is a function, and we need to treat the data as if we had plotted
136	c it on a graph as points ( i.e it is a ratio not a histogram) so:
137	c 1. determine the values of Q, the midpoints of the X bin boundaries.
138	c 2. transform Y and E to Y' = F(Y) and E' = E*(dF(Y)/dY).
139	c this error calc. assumes the errors are uncorrelated, which is not
140	c actually true e.g. if F(Y)=Y**2 E'= 2.E.F(Y)
141	c wheras one usually takes E'= SQRT(2).E.F(Y)
142	c 3. transform X to X' = f(X).
143	c (do not multiply Y by dX !)
144	c
145	c the rules for histograms are :
146	c 1. multiply Y and E by the Jacobian J(x,x') = dx/dx' which converts
147	c the data from bins in x to bins in x'.
148	c 2. transform Y and E to Y' = F(Y) and E' = E*(dF(Y)/dY).
149	c 3. transform X to X' = f(X).
150	c (do not multiply Y by dX !)
151	c
152	c if(option.eq.4.or.option.eq.104.or.option.eq.6.or.option.eq.106)then
153	cc==== REVERSE the order of the data first
154	c do i = 1,npt(n)
155	c x(n2 + 2 - i) = x(offset_n + i)
156	c y(n2 + 1 - i)=(y(offset_n+i) - bkgd)
157	c e(n2 + 1 - i) = e(offset_n+i)
158	c end do
159	c x(offset_i_graph) = x(offset_n+npt(n)+1)
160	c else
161
162	* FIRST TAKE OFF THE BACKGROUND
163	if (subtract_wkspc.eq.0) then
164	* FOR SINGLE VALUE
165	do i = 1,npt(n)
166	x(offset_i_graph+i)=x(offset_n+i)
167	y(offset_i_graph+i)=y(offset_n+i) - bkgd
168	e(offset_i_graph+i)=e(offset_n+i)
169	end do
170	x(offset_i_graph+npt(n)+1) = x(offset_n+npt(n)+1)
171	else if (subtract_wkspc.eq.1) then
172	* FOR BACKGROUND IN WORKSPACE
173	if (npt(n).ne.npt(j_graph)) then
174	call error('ERROR: Data & Background workspaces have
175	& different numbers of points. Rebin. ')
176	goto 102
177	end if
178	do i = 1,npt(n)
179	if (x(offset_n+i).ne.x(offset_j_graph+i)) then
180	call error('ERROR: Data & Background workspaces have
181	& different X scales. Rebin. ')
182	goto 102
183	end if
184	x(offset_i_graph+i)=x(offset_n+i)
185	y(offset_i_graph+i)=y(offset_n+i) - y(offset_j_graph+i)
186	e(offset_i_graph+i)=e(offset_n+i)
187	end do
188	x(offset_i_graph+npt(n)+1) = x(offset_n+npt(n)+1)
189	end if
190
191	c end if
192	c
193	c==== Guinier Plot ln(I) vs Q**2
194	if ( option .eq. 1 ) then
195	do i = n1,n2
196	if (y(i) .gt. 0.0) then
197	e(i) = e(i)/y(i)
198	y(i) = log( y(i))
199	else
200	y(i) = 0.0
201	e(i) = 0.0
202	end if
203	x(i) = x(i)**2
204	end do
205	x(n2+1)=x(n2+1)**2
206	xcaption(i_graph) = 'Q\u2\d (\A\u-2\d)'
207	ycaption(i_graph) = 'log\de\u(I)'
208
209	c==== Guinier Plot (rods) ln(IQ) vs Q**2
210
211	else if (option .eq. 2) then
212	c
213	do i=n1,n2
214	if (y(i) .gt. 0.0) then
215
216	if(hist)then
217	Q = (x(i)+x(i+1))*0.5
218	else
219	Q = x(i)
220	end if
221
222	e(i) = e(i)/y(i)
223	y(i) = log( y(i)*Q)
224	else
225	y(i) = 0.0
226	e(i) = 0.0
227	end if
228	x(i) = x(i)**2
229	end do
230	x(n2+1)=x(n2+1)**2
231	xcaption(i_graph) = 'Q\u2\d (\A\u-2\d)'
232	ycaption(i_graph) = 'log\de\u(IQ)'
233
234	c==== Guinier Plot (sheets) ln(I.Q2) vs Q2
235
236	else if (option .eq. 3) then
237	c
238	do i=n1,n2
239	if (y(i) .gt. 0.0) then
240
241	if(hist)then
242	Q = (x(i)+x(i+1))*0.5
243	else
244	Q = x(i)
245	end if
246
247	e(i) = e(i)/y(i)
248	y(i) = log( y(i)Q*2)
249	else
250	y(i) = 0.0
251	e(i) = 0.0
252	end if
253	x(i) = x(i)**2
254	end do
255	x(n2+1)=x(n2+1)**2
256	xcaption(i_graph) = 'Q\u2\d (\A\u-2\d)'
257	ycaption(i_graph) = 'log\de\u(I.Q\u2\d)'
258
259	c==== Zimm Plot 1/I vs Q**2
260
261	else if (option .eq. 4) then
262	c
263	do i=n1,n2
264	if(y(i).gt.0.0)then
265	y(i) = 1.0 / y(i)
266	e(i) = e(i) * y(i)**2
267	else
268	e(i)=0.0
269	y(i)=0.0
270	end if
271	x(i) = x(i)**2
272	end do
273	x(n2+1)=x(n2+1)**2
274	xcaption(i_graph) = 'Q\u2\d (\A\u-2\d)'
275	ycaption(i_graph) = 'I\u-1\d'
276
277	c==== Kratky Plot I.Q2 vs Q (changed from Q2, RKH 15/4/92)
278
279	else if (option .eq. 5) then
280	c
281	do i=n1,n2 ! multiply y by Q**2
282
283	if(hist)then
284	Q = (x(i)+x(i+1))*0.5
285	else
286	Q = x(i)
287	end if
288
289	y(i) = y(i)Q*2
290	e(i) = e(i)Q*2
291	end do
292	xcaption(i_graph) = 'Q (\A\u-1\d)'
293	ycaption(i_graph) = 'I.Q\u2\d'
294	c
295	c==== Debye-Bueche Plot 1/sqrt(I) vs Q**2
296	c
297	else if (option .eq. 6) then
298	c
299	do i = n1,n2
300	if (y(i) .gt. 0.0) then
301	y(i) = 1 / sqrt (y(i))
302	e(i) = e(i)0.5 y(i)**3
303	else
304	y(i) = 0.0
305	e(i) = 0.0
306	end if
307	x(i) = x(i)**2
308	end do
309	x(n2+1)=x(n2+1)**2
310	xcaption(i_graph) = 'Q\u2\d (\A\u-2\d)'
311	ycaption(i_graph) = 'I\u-\(261)\d'
312	c
313	c==== LOG - LOG plot
314	else if (option .eq. 7)then
315	ii=n2
316	if(hist)ii=n2+1
317	do i=n1,ii
318	if (x(i) .gt. 0) then
319	x(i)= log(x(i))
320	else
321	call error('ERROR: log of negative or zero in x array')
322	get_special = status_not_ok
323	call replace_prompt
324	return
325	end if
326	end do
327	c
328	do i= n1,n2
329	if (y(i) .gt. 0) then
330	e(i) = e(i) / y(i)
331	y(i) = log( y(i))
332	else
333	y(i) = 0
334	e(i) = 0
335	endif
336	end do
337	xcaption(i_graph) = 'log\de\u(Q)'
338	ycaption(i_graph) = 'log\de\u(I)'
339	c
340	c==== Porod Plot I.Q**4 vs Q
341	c
342	else if (option .eq. 8) then
343	c
344	do i = n1,n2
345
346	if(hist)then
347	Q = (x(i)+x(i+1))*0.5
348	else
349	Q = x(i)
350	end if
351
352	y(i) = y(i)Q*4
353	e(i) = e(i)Q*4
354	end do
355	c
356	xcaption(i_graph) = 'Q (\A\u-1\d)'
357	ycaption(i_graph) = 'I.Q\u4\d'
358	c
359	c=== general function (QA)(IB)xLn(QCxID) v similar',
360	else if( option .eq. 9 .or. option .eq. 91) then
361	if (option.eq. 9)then
362	status=smg$save_physical_screen(pbid,save_vdid,4)
363	write(6,1001)
364	1001 format(1x,'Special function plots',/,
365	* ' Q and I are transformed by equations of the type:',/,
366	* ' Q = Q^A * I^B * LOGe( Q^C * I^D * E)',/,
367	* ' I = Q^F * I^G * LOGe( Q^H * I^I * J)',/,
368	* ' (Note: LOGe(0.0) is taken as 1.0 )',/,
369	* ' Please enter A, B, C, D & E >>> ',$)
370	read(5,*)SP_E,SP_F,SP_G,SP_H,SP_CX
371	write(6,1002)
372	1002 format(/,' Please enter F, G, H, I & J >>> ',$)
373	read(5,*)SP_A,SP_B,SP_C,SP_D,SP_CY
374	c23456789012345678901234567890123456789012345678901234567890123456789012345678901234567890
375	write(6,1003)SP_E,SP_F,SP_G,SP_H,SP_CX,SP_A,SP_B,SP_C,SP_D,SP_CY
376	1003 format(1x,' CHECK: Q = Q*',f6.3,' I**',f6.3,
377	& ' * LOGe{ Q*',f6.3,' I*',f6.3,' ',f6.3,' }',/,
378	& 1x,' I = Q*',f6.3,' I**',f6.3,
379	& ' * LOGe{ Q*',f6.3,' I*',f6.3,' ',f6.3,' }',/,
380	& 2x,' OK [Answer 1=YES, 0=NO]? ')
381	read(5,*)i
382	if(i.ne.1)then
383	get_special = status_not_ok
384	status=smg$restore_physical_screen(pbid,save_vdid)
385	call replace_prompt
386	return
387	end if
388	c
389	else
390	call remark(' using previous special function')
391	end if
392	c==== store the original mid point x values and transform the old ones
393	c==== function PW(a,b) does a**b, function Slog does LOG( ) both
394	c==== checking for errors where possible.
395	offset = len_streams * (used_streams - (i_graph-1))
396	c==== new x(i) values are transform of old x(i) and old y(i)
397	c==== except last point of histogram where use last old x(n2+1), y(n2)
398	c==== sorted out a muddle in the offsets here RKH 29/10/92
399	do i = n1 , n2
400	Q=x(i)
401	x(i) = PW(y(i),SP_F) * PW(Q,SP_E) *
402	> Slog( PW(y(i),SP_H) * PW(Q,SP_G) * SP_CX )
403	end do
404	if(hist)x(n2+1) = PW(y(n2),SP_F) * PW( x(n2+1),SP_E ) *
405	> Slog( PW(y(n2),SP_H) * PW( x(n2+1),SP_G )* SP_CX ) ! fix up the last point
406	c
407	c==== new y(i) values are transform of old midpoint ( x(i), x(i+1) )
408	c==== if histogram,[ or just x(i) if not ] with old y(i)
409	ii=0
410	do i= n1,n2
411	ii=ii+1
412	c==== care to point at the OLD x values !
413	if(hist)then
414	Q = (x(offset_n + ii)+x(offset_n + ii+1))*0.5
415	else
416	Q = x(offset_n +ii)
417	end if
418
419	XCYD = PW( y(i),SP_D ) * PW( Q,SP_C ) * SP_CY
420	RLG = Slog( XCYD )
421	e(i) = e(i)* PW(Q,SP_A) * ( DERV(y(i),SP_B) * RLG +
422	> PW(y(i),SP_B) * PW(XCYD,-1.0) * PW( Q,SP_C ) *
423	> SP_CY * DERV(y(i),SP_D) )
424
425	y(i) = PW( y(i),SP_B ) * PW( Q,SP_A ) * RLG
426	end do
427
428	status=smg$restore_physical_screen(pbid,save_vdid)
429	xcaption(i_graph) = 'Special'
430	ycaption(i_graph) = 'Special'
431
432	c======================= SECTION FOR HISTOGRAM DATA ===========
433	c=====do not multiply by dX first!
434	c note Q is given by Q(i) = (x(i)+x(i+1))/2.
435	c QJ is the Jacobian = dx/dx'. For x=Q and x'=Q**2, QJ=1/2Q.
436	c==== Guinier Plot
437	else if (option .eq. 101) then
438	c
439	do i = n1,n2
440	QJ = 1/(x(i) + x(i+1))
441	y(i) = y(i)*QJ
442	e(i) = e(i)*QJ
443	if (y(i) .gt. 0.0) then
444	e(i) = e(i)/y(i)
445	y(i) = log( y(i))
446	else
447	y(i) = 0.0
448	e(i) = 0.0
449	end if
450	end do
451	c
452	do i=n1,n2+1 ! convert x to Q**2
453	x(i) = x(i) ** 2
454	end do
455	xcaption(i_graph) = 'Q\u2\d (\A\u-2\d)'
456	ycaption(i_graph) = 'log\de\u(I)'
457
458	c==== Guinier Plot (rods)
459
460	else if (option .eq.102) then
461	c
462	do i = n1,n2
463	QJ = 1/(x(i) + x(i+1))
464	y(i) = y(i)*QJ
465	e(i) = e(i)*QJ
466	if (y(i) .gt. 0.0) then
467	e(i) = e(i)/y(i)
468	y(i) = log( y(i)/(2*QJ))
469	else
470	y(i) = 0.0
471	e(i) = 0.0
472	end if
473	end do
474	c
475	do i=n1,n2+1 ! convert x to Q**2
476	x(i) = x(i) ** 2
477	end do
478	xcaption(i_graph) = 'Q\u2\d (\A\u-2\d)'
479	ycaption(i_graph) = 'log\de\u(I.Q)'
480
481	c==== Guinier Plot (sheets)
482
483	else if (option .eq. 103) then
484	c
485	do i = n1,n2
486	QJ = 1/(x(i) + x(i+1))
487	y(i) = y(i)*QJ
488	e(i) = e(i)*QJ
489	if (y(i) .gt. 0.0) then
490	e(i) = e(i)/y(i)
491	y(i) = log( y(i)/(4QJQJ))
492	else
493	y(i) = 0.0
494	e(i) = 0.0
495	end if
496	end do
497	c
498	do i=n1,n2+1 ! convert x to Q**2
499	x(i) = x(i) ** 2
500	end do
501	xcaption(i_graph) = 'Q\u2\d (\A\u-2\d)'
502	ycaption(i_graph) = 'log\de\u(I.Q\u2\d)'
503
504	c==== Zimm Plot 1/I vs Q**2
505
506	else if (option .eq. 104) then
507	c
508	do i=n1,n2
509	if(y(i).gt.0.0)then
510	QJ = 1/(x(i) + x(i+1))
511	y(i) = y(i)*QJ
512	e(i) = e(i)*QJ
513	y(i) = 1/y(i)
514	e(i) = e(i)y(i)*2
515	else
516	e(i)=0.0
517	y(i)=0.0
518	end if
519	end do
520	c
521	do i=n1,n2+1
522	x(i)=x(i)**2
523	end do
524	xcaption(i_graph) = 'Q\u2\d (\A\u-2\d)'
525	ycaption(i_graph) = 'I\u-1\d'
526
527	c==== Kratky Plot I.Q2 vs Q (changed from Q2 RKH 15/4/92)
528
529	else if (option .eq. 105) then
530	c
531	do i=n1,n2 ! multiply y by Q**2
532	c QJ = 1/(x(i) + x(i+1))
533	c y(i) = y(i)*QJ
534	c e(i) = e(i)*QJ
535	c y(i) = y(i)/(4QJQJ)
536	c e(i) = e(i)/(4QJQJ)
537	QJ = 0.5*(x(i) + x(i+1))
538	y(i) = y(i)QJQJ
539	e(i) = e(i)QJQJ
540	end do
541	c
542	c do i=n1,n2+1
543	c x(i)=x(i)**2
544	c end do
545	xcaption(i_graph) = 'Q (\A\u-1\d)'
546	ycaption(i_graph) = 'I.Q\u2\d'
547	c
548	c==== Debye-Bueche Plot 1/sqrt(I) vs Q**2
549	c
550	else if (option .eq. 106) then
551	c
552	do i=n1,n2
553	QJ = 1/(x(i) + x(i+1))
554	y(i) = y(i)*QJ
555	e(i) = e(i)*QJ
556	if(y(i).gt.0.0)then
557	y(i) = 1/sqrt(y(i))
558	e(i) = e(i)0.5y(i)**3
559	else
560	e(i)=0.0
561	y(i)=0.0
562	end if
563	end do
564	c
565	do i=n1,n2+1
566	x(i)=x(i)**2
567	end do
568	xcaption(i_graph) = 'Q\u2\d (\A\u-2\d)'
569	ycaption(i_graph) = 'I.Q\u2\d'
570	c
571	c==== LOG - LOG plot
572	else if (option .eq. 107)then
573	c
574	do i= n1,n2
575	Q = (x(i) + x(i+1))/2
576	y(i) = y(i)/Q
577	e(i) = e(i)/Q
578	if(y(i).gt.0)then
579	e(i) = e(i)/y(i)
580	y(i) = log( y(i))
581	else
582	y(i)=0
583	e(i)=0
584	endif
585	end do
586	c
587	do i=n1,n2+1
588	if(x(i).gt.0)then
589	x(i)= log(x(i))
590	else
591	call error('ERROR: log of negative or zero in x array')
592	get_special = status_not_ok
593	call replace_prompt
594	return
595	end if
596	end do
597	xcaption(i_graph) = 'log\de\u(Q)'
598	ycaption(i_graph) = 'log\de\u(I)'
599
600	c==== Porod Plot I.Q4 vs Q2
601	c
602	else if (option .eq. 108) then
603	c
604	do i = n1,n2
605	Q = (x(i) + x(i+1))/2
606	y(i) = y(i)Q*4
607	e(i) = e(i)Q*4
608	end do
609	c
610	do i=n1,n2+1
611	x(i)=x(i)**2
612	end do
613	xcaption(i_graph) = 'Q\u2\d (\A\u-2\d)'
614	ycaption(i_graph) = 'I.Q\u4\d'
615	c
616	c=== general function (QA)(IB)xLn(QCxID) v similar',
617	else if( option .eq. 109) then
618	status=smg$save_physical_screen(pbid,save_vdid,4)
619	write(6,1001)
620	read(5,*)A,B,C,D,CY
621	write(6,1002)
622	read(5,*)EE,F,G,H,CX
623	write(6,1003)A,B,C,D,CY,EE,F,G,H,CX
624	read(5,*)I
625	if(I.ne.1)then
626	get_special = status_not_ok
627	call replace_prompt
628	return
629	end if
630	c==== store the original x values and transform the old ones
631	c==== function PW(a,b) does a**b, function Slog does LOG( ) both
632	c==== checking for errors where possible.
633	offset = len_streams * (used_streams - (i_graph-1))
634	do i = n1 , n2
635	x(i+offset) =0.5*( x(i+1) + x(i))
636	yy= 0.5*(y(i+1)+y(i))
637	x(i) = PW(yy,F) * PW(x(i),EE) *
638	> Slog( PW(yy,H) * PW(x(i),G) * CX )
639	end do
640	x(n2+1) = PW(yy,F) * PW(x(n2+1),EE) * Slog( PW(yy,H)
641	> * PW(x(n2+1),G)* CX ) ! fix up the last point
642	c
643	c==== now rescale the old Y and E before transforming them
644	do i=n1,n2
645	y(i) = y(i) / (x(i+1)-x(i))
646	e(i) = e(i) / (x(i+1)-x(i))
647	XCYD = PW(y(i),D) * PW(x(i+offset),C) * CY
648	RLG = Slog( XCYD )
649	e(i) = e(i)* PW(x(i+offset),A) * ( DERV(y(i),B) * RLG +
650	> PW(y(i),B) * PW(XCYD,-1.0) * PW( x(i+offset),C ) *
651	> CY * DERV(y(i),D) )
652
653	y(i) = PW(y(i),B) * PW(x(i+offset),A) * RLG
654	end do
655
656	status=smg$restore_physical_screen(pbid,save_vdid)
657	xcaption(i_graph) = 'Special'
658	ycaption(i_graph) = 'Special'
659
660	end if
661
662	xmn(i_graph) = x(n1)
663	xmx(i_graph) = x(n2)
664	if(hist)xmx(i_graph) = x(n2+1)
665	ymn(i_graph) = y(n1)
666	ymx(i_graph) = y(n1)
667	do i = n1,n2
668	ymn(i_graph) = min (ymn(i_graph), y(i))
669	ymx(i_graph) = max (ymx(i_graph), y(i))
670	end do
671
672	c===============================================================================
673	c copy header into workspace i_graph
674	c
675	long_title(i_graph) = long_title(n)
676	run_number(i_graph) = run_number(n)
677	dat_file(i_graph) = dat_file(n)
678	run_user(i_graph) = run_user(n)
679	start_time(i_graph) = start_time(n)
680	bin_parameter(i_graph) = bin_parameter(n)
681	len_pdfn(i_graph) = len_pdfn(n)
682
683	c===============================================================================
684	c Fill workspace
685
686	write (text, '(i5)') n
687	work(i_graph).label_1 = 'WORKSPACE : '//text
688	if (work(n).lambda_max .gt. 0) then
689	write (text, '(f5.2, a, f5.2, a)')
690	> lambda_min, ' > ', lambda_max, ' Ang'
691	work(i_graph).label_2 = 'LAMBDAS : '//text
692	else
693	work(i_graph).label_2 = 'All wavelengths'
694	end if
695	write (text, '(f5.0, a, f5.0, a)') phi_min,' > ',phi_max, ' deg'
696	work(i_graph).label_3 = 'PHI : '//text
697	xcode(i_graph) = 7
698	ycode(i_graph) = xcode(i_graph)
699
700	work(i_graph).nsp_first = work(n).nsp_first
701	work(i_graph).nsp_last = work(n).nsp_last
702	work(i_graph).lambda_min = work(n).lambda_min
703	work(i_graph).lambda_max = work(n).lambda_max
704	work(i_graph).t_min = work(n).t_min
705	work(i_graph).t_max = work(n).t_max
706	work(i_graph).r_min = work(n).r_min
707	work(i_graph).r_max = work(n).r_max
708	work(i_graph).sample_file = work(n).sample_file
709	work(i_graph).can_file = work(n).can_file
710	work(i_graph).norm_file = work(n).norm_file
711	work(i_graph).back_file = work(n).back_file
712
713	get_special = status_ok
714
715	goto 101
716
717	100 call error ('ERROR: Invalid entry')
718
719	102 get_special = status_not_ok
720
721	101 call replace_prompt
722
723	return
724
725	end
726	c
727	c===============================================================================
728	c COLETTE - LOQ data analysis development program - R.Heenan 3/9/87
729	c
730	c=== this function does A**B, trying to avoid error returns
731	c=== it doe NOT try to anticipate overflows however.
732	c
733	function pw(a,b)
734	c
735	if(abs(a).lt.1.0e-30)then ! 0n or 0-n are set zero
736	pw=0.0
737	c
738	else if (abs(b).lt.1.0e-08)then ! a**0 is set to 1.0
739	pw=1.0
740	c
741	else if(a.lt.0.0)then ! test for (-a)**b
742	c
743	if( abs( b - nint(b) ).lt.1.e-08 ) then
744	if(nint(b).eq.1)then
745	pw = a ! (-a)**1
746	else if(nint(b).eq.-1)then
747	pw = 1./a ! (-a)**(-1)
748	else
749	if(mod(nint(b),2).eq.0)then
750	pw=(-a)**b
751	else
752	pw=-(-a)**b
753	end if
754	end if
755	return
756	else
757	write(6,1005)
758	1005 format(1x,'ERROR: negative number to non-integral power')
759	pw=0.0
760	return
761	end if
762	return
763	c
764	else
765	pw=a**b ! the "normal" route
766	c
767	end if
768	return
769	c
770	end
771	c
772	c
773	c===============================================================================
774	c COLETTE - LOQ data analysis development program - R.Heenan 3/9/87
775	c
776	c=== this function does the derivative of A**B, with respect to A
777	c=== trying to avoid error returns
778	c=== it doe NOT try to anticipate overflows however.
779	c
780	function derv(a,b)
781	c
782	if( abs(a).lt.1.0e-30 ) then
783	derv = 0.0
784	c
785	else if ( abs(b) .lt.1.0e-08 ) then
786	derv = 0.0
787	c
788	else
789	derv = b* PW( a, b-1.0 )
790	end if
791	c
792	return
793	end
794	c
795	c===============================================================================
796	c COLETTE - LOQ data analysis development program - R.Heenan 3/9/87
797	c
798	c=== this function does LOGe(a), trying to avoid error returns
799	c
800	function Slog(a)
801	if(a.gt.0.0)then
802	slog = log(a)
803	else
804	slog = 1.0
805	end if
806	return
807	end
808	c

Download in other formats:

Original Format