OpenCMISS-Iron Internal API Documentation
coordinate_routines.f90
Go to the documentation of this file.
1 
43 
45 MODULE coordinate_routines
46 
47  USE base_routines
48  USE constants
49  USE input_output
50  USE iso_varying_string
51  USE kinds
52  USE maths
53  USE strings
54  USE types
55 
56 #include "macros.h"
57 
58  IMPLICIT NONE
59 
60 !!TODO: Should all of the get/set/create/destroy routines be accessed via pointers???
61 
62  PRIVATE
63 
64  !Module parameters
65 
70  INTEGER(INTG), PARAMETER :: coordinate_rectangular_cartesian_type=1
71  INTEGER(INTG), PARAMETER :: coordinate_cylindrical_polar_type=2
72  INTEGER(INTG), PARAMETER :: coordinate_spherical_polar_type=3
73  INTEGER(INTG), PARAMETER :: coordinate_prolate_spheroidal_type=4
74  INTEGER(INTG), PARAMETER :: coordinate_oblate_spheroidal_type=5
76 
81  INTEGER(INTG), PARAMETER :: coordinate_no_radial_interpolation_type=0
82  INTEGER(INTG), PARAMETER :: coordinate_radial_interpolation_type=1
83  INTEGER(INTG), PARAMETER :: coordinate_radial_squared_interpolation_type=2
84  INTEGER(INTG), PARAMETER :: coordinate_radial_cubed_interpolation_type=3
86 
91  INTEGER(INTG), PARAMETER :: coordinate_jacobian_no_type=0
92  INTEGER(INTG), PARAMETER :: coordinate_jacobian_line_type=1
93  INTEGER(INTG), PARAMETER :: coordinate_jacobian_area_type=2
94  INTEGER(INTG), PARAMETER :: coordinate_jacobian_volume_type=3
96 
97  !Module types
98 
99  TYPE coordinate_systems_type
100  INTEGER(INTG) :: number_of_coordinate_systems
101  TYPE(coordinate_system_ptr_type), POINTER :: coordinate_systems(:)
102  END TYPE coordinate_systems_type
103 
104  !Module variables
105 
106  CHARACTER(LEN=21) :: coordinate_system_type_string(5) = &
107  & (/ "Rectangular Cartesian",&
108  & "Cylindrical Polar ", &
109  & "Spherical Polar ", &
110  & "Prolate Spheroidal ", &
111  & "Oblate Spheroidal " /)
112 
113  TYPE(coordinate_systems_type) :: coordinate_systems
114 
115  !Interfaces
116 
119  INTERFACE coordinate_convert_from_rc
120  MODULE PROCEDURE coordinate_convert_from_rc_dp
121  MODULE PROCEDURE coordinate_convert_from_rc_sp
122  END INTERFACE !COORDINATE_CONVERT_FROM_RC
123 
126  INTERFACE coordinate_convert_to_rc
127  MODULE PROCEDURE coordinate_convert_to_rc_dp
128  MODULE PROCEDURE coordinate_convert_to_rc_sp
129  END INTERFACE !COORDINATE_CONVERT_TO_RC
130 
133  INTERFACE coordinate_delta_calculate
134  MODULE PROCEDURE coordinate_delta_calculate_dp
135  !MODULE PROCEDURE COORDINATE_DELTA_CALCULATE_SP
136  END INTERFACE !COORDINATE_DELTA_CALCULATE
137 
139  INTERFACE dxz
140  MODULE PROCEDURE dxz_dp
141  !MODULE PROCEDURE DXZ_SP
142  END INTERFACE !Dxz
143 
144  !!TODO:: CHANGE NAME TO SOMETHING MORE MEANINGFULL?
145  INTERFACE d2zx
146  MODULE PROCEDURE d2zx_dp
147  !MODULE PROCEDURE D2ZX_SP
148  END INTERFACE !D2ZX
149 
150  !!TODO:: CHANGE NAME TO SOMETHING MORE MEANINGFULL?
151  INTERFACE dzx
152  MODULE PROCEDURE dzx_dp
153  !MODULE PROCEDURE DZX_SP
154  END INTERFACE !DZX
155 
156  INTERFACE coordinate_derivative_convert_to_rc
157  MODULE PROCEDURE coordinate_derivative_convert_to_rc_dp
158  MODULE PROCEDURE coordinate_derivative_convert_to_rc_sp
159  END INTERFACE !COORDINATE_DERIVATIVE_CONVERT_TO_RC
160 
161  PUBLIC coordinate_rectangular_cartesian_type,coordinate_cylindrical_polar_type,coordinate_spherical_polar_type, &
162  & coordinate_prolate_spheroidal_type,coordinate_oblate_spheroidal_type
163 
164  PUBLIC coordinate_no_radial_interpolation_type,coordinate_radial_interpolation_type, &
165  & coordinate_radial_squared_interpolation_type,coordinate_radial_cubed_interpolation_type
166 
167  PUBLIC coordinate_jacobian_no_type,coordinate_jacobian_line_type,coordinate_jacobian_area_type,coordinate_jacobian_volume_type
168 
169  PUBLIC coordinate_system_type_string
170 
171  PUBLIC coordinate_convert_from_rc,coordinate_convert_to_rc,coordinate_delta_calculate,coordinate_derivative_norm, &
172  & coordinate_interpolation_adjust,coordinate_interpolation_parameters_adjust,coordinate_metrics_calculate
173 
174  PUBLIC coordinate_system_dimension_get,coordinate_system_dimension_set
175 
176  PUBLIC coordinate_system_focus_get,coordinate_system_focus_set
177 
178  PUBLIC coordinates_radialinterpolationtypeget,coordinates_radialinterpolationtypeset
179 
180  PUBLIC coordinate_system_type_get,coordinate_system_type_set
181 
182  PUBLIC coordinate_system_origin_get,coordinate_system_origin_set
183 
184  PUBLIC coordinate_system_orientation_get,coordinate_system_orientation_set
185 
186  PUBLIC coordinate_system_create_start,coordinate_system_create_finish
187 
188  PUBLIC coordinate_system_destroy
189 
190  PUBLIC coordinate_derivative_convert_to_rc
191 
192  PUBLIC coordinate_system_user_number_find
193 
194  PUBLIC coordinate_systems_initialise,coordinate_systems_finalise
195 
196  PUBLIC coordinates_materialsystemcalculate
197 
198 CONTAINS
199 
200  !
201  !================================================================================================================================
202  !
203 
205  FUNCTION coordinate_convert_from_rc_dp(COORDINATE_SYSTEM,Z,ERR,ERROR)
206 
207  !Argument variables
208  TYPE(coordinate_system_type), INTENT(IN) :: COORDINATE_SYSTEM
209  REAL(DP), INTENT(IN) :: Z(:)
210  INTEGER(INTG), INTENT(OUT) :: ERR
211  TYPE(varying_string), INTENT(OUT) :: ERROR
212  !Function variable
213  REAL(DP) :: COORDINATE_CONVERT_FROM_RC_DP(size(z,1))
214  !Local variables
215  REAL(DP) :: A1,A2,A3,A4,A5,A6,A7,A8,A9,FOCUS
216 
217  enters("COORDINATE_CONVERT_FROM_RC_DP",err,error,*999)
218 
219  coordinate_convert_from_rc_dp=0.0_dp
220 
221  IF(SIZE(z,1)<coordinate_system%NUMBER_OF_DIMENSIONS) &
222  & CALL flagerror("Size of Z is less than the number of dimensions.",err,error,*999)
223 
224  SELECT CASE(coordinate_system%TYPE)
225  CASE(coordinate_rectangular_cartesian_type)
226  coordinate_convert_from_rc_dp(1:coordinate_system%NUMBER_OF_DIMENSIONS)=z(1:coordinate_system%NUMBER_OF_DIMENSIONS)
227  CASE(coordinate_cylindrical_polar_type)
228  SELECT CASE(coordinate_system%NUMBER_OF_DIMENSIONS)
229  CASE(2)
230  coordinate_convert_from_rc_dp(1)=sqrt(z(1)**2+z(2)**2)
231  coordinate_convert_from_rc_dp(2)=atan2(z(1),z(2))
232  CASE(3)
233  coordinate_convert_from_rc_dp(1)=sqrt(z(1)**2+z(2)**2)
234  coordinate_convert_from_rc_dp(2)=atan2(z(1),z(2))
235  coordinate_convert_from_rc_dp(3)=z(3)
236  CASE DEFAULT
237  CALL flagerror("Invalid number of coordinates.",err,error,*999)
238  END SELECT
239  IF(coordinate_convert_from_rc_dp(2)<0.0_dp) &
240  & coordinate_convert_from_rc_dp(2)=coordinate_convert_from_rc_dp(2)+2.0_dp*pi !reference coordinate 0->2*pi
241  CASE(coordinate_spherical_polar_type)
242  IF(coordinate_system%NUMBER_OF_DIMENSIONS==3) THEN
243  coordinate_convert_from_rc_dp(1)=sqrt(z(1)**2+z(2)**2+z(3)**2)
244  IF(abs(z(1))>=zero_tolerance.OR.abs(z(2))>=zero_tolerance) THEN
245  coordinate_convert_from_rc_dp(2)=atan2(z(2),z(1))
246  ELSE
247  coordinate_convert_from_rc_dp(2)=0.0_dp
248  ENDIF
249  a1=sqrt(z(1)**2+z(2)**2)
250  IF(abs(z(3))>=zero_tolerance.OR.abs(a1)>=zero_tolerance) THEN
251  coordinate_convert_from_rc_dp(3)=atan2(z(3),a1)
252  ELSE
253  coordinate_convert_from_rc_dp(3)=0.0_dp
254  ENDIF
255  ELSE
256  CALL flagerror("Invalid number of coordinates.",err,error,*999)
257  ENDIF
258  CASE(coordinate_prolate_spheroidal_type)
259  focus=coordinate_system%FOCUS
260  IF(coordinate_system%NUMBER_OF_DIMENSIONS==3) THEN
261  a1=z(1)**2+z(2)**2+z(3)**2-focus**2
262  a2=sqrt(a1**2+4.0_dp*(focus**2)*(z(2)**2+z(3)**2))
263  a3=2.0_dp*focus**2
264  a4=max((a2+a1)/a3,0.0_dp)
265  a5=max((a2-a1)/a3,0.0_dp)
266  a6=sqrt(a4)
267  a7=min(sqrt(a5),1.0_dp)
268  IF(abs(a7)<=1.0_dp) THEN
269  a8=asin(a7)
270  ELSE
271  a8=0.0_dp
272  CALL flag_warning("Put A8=0 since ABS(A8)>1.",err,error,*999)
273  ENDIF
274  IF((abs(z(3))<zero_tolerance).OR.(abs(a6)<zero_tolerance).OR.(abs(a7)<zero_tolerance)) THEN
275  a9=0.0_dp
276  ELSE
277  IF(abs(a6*a7)>0.0_dp) THEN
278  a9=z(3)/(focus*a6*a7)
279  ELSE
280  a9=0.0_dp
281  CALL flag_warning("Put A9=0 since A6*A7=0.",err,error,*999)
282  ENDIF
283  IF(a9>=1.0_dp) THEN
284  a9=pi/2.0_dp
285  ELSE IF(a9<=-1.0_dp) THEN
286  a9=-pi/2.0_dp
287  ELSE
288  a9=asin(a9)
289  ENDIF
290  ENDIF
291  coordinate_convert_from_rc_dp(1)=log(a6+sqrt(a4+1.0_dp))
292  IF(abs(z(1))>=zero_tolerance) THEN
293  coordinate_convert_from_rc_dp(2)=a8
294  ELSE
295  coordinate_convert_from_rc_dp(2)=pi-a8
296  ENDIF
297  IF(abs(z(2))>zero_tolerance) THEN
298  coordinate_convert_from_rc_dp(3)=mod(a9+2.0_dp*pi,2.0_dp*pi)
299  ELSE
300  coordinate_convert_from_rc_dp(3)=pi-a9
301  ENDIF
302  ELSE
303  CALL flagerror("Invalid number of coordinates.",err,error,*999)
304  ENDIF
305  CASE(coordinate_oblate_spheroidal_type)
306  CALL flagerror("Not implemented.",err,error,*999)
307  CASE DEFAULT
308  CALL flagerror("Invalid coordinate type.",err,error,*999)
309  END SELECT
310 
311  exits("COORDINATE_CONVERT_FROM_RC_DP")
312  RETURN
313 999 errorsexits("COORDINATE_CONVERT_FROM_RC_DP",err,error)
314  RETURN
315  END FUNCTION coordinate_convert_from_rc_dp
316 
317  !
318  !================================================================================================================================
319  !
320 
324  FUNCTION coordinate_convert_from_rc_sp(COORDINATE_SYSTEM,Z,ERR,ERROR)
325 
326  !Argument variables
327  TYPE(coordinate_system_type), INTENT(IN) :: COORDINATE_SYSTEM
328  REAL(SP), INTENT(IN) :: Z(:)
329  INTEGER(INTG), INTENT(OUT) :: ERR
330  TYPE(varying_string), INTENT(OUT) :: ERROR
331  !Function variable
332  REAL(SP) :: COORDINATE_CONVERT_FROM_RC_SP(size(z,1))
333  !Local variables
334  REAL(SP) :: A1,A2,A3,A4,A5,A6,A7,A8,A9,FOCUS
335 
336  enters("COORDINATE_CONVERT_FROM_RC_SP",err,error,*999)
337 
338  coordinate_convert_from_rc_sp=0.0_sp
339 
340  IF(SIZE(z,1)<coordinate_system%NUMBER_OF_DIMENSIONS) &
341  & CALL flagerror("Size of Z is less than the number of dimensions.",err,error,*999)
342 
343  SELECT CASE(coordinate_system%TYPE)
344  CASE(coordinate_rectangular_cartesian_type)
345  coordinate_convert_from_rc_sp(1:coordinate_system%NUMBER_OF_DIMENSIONS)=z(1:coordinate_system%NUMBER_OF_DIMENSIONS)
346  CASE(coordinate_cylindrical_polar_type)
347  SELECT CASE(coordinate_system%NUMBER_OF_DIMENSIONS)
348  CASE(2)
349  coordinate_convert_from_rc_sp(1)=sqrt(z(1)**2+z(2)**2)
350  coordinate_convert_from_rc_sp(2)=atan2(z(1),z(2))
351  CASE(3)
352  coordinate_convert_from_rc_sp(1)=sqrt(z(1)**2+z(2)**2)
353  coordinate_convert_from_rc_sp(2)=atan2(z(1),z(2))
354  coordinate_convert_from_rc_sp(3)=z(3)
355  CASE DEFAULT
356  CALL flagerror("Invalid number of coordinates.",err,error,*999)
357  END SELECT
358  IF(coordinate_convert_from_rc_sp(2)<0.0_sp) &
359  & coordinate_convert_from_rc_sp(2)=coordinate_convert_from_rc_sp(2)+2.0_sp*REAL(PI,SP) !reference coordinate 0->2*pi
360  CASE(coordinate_spherical_polar_type)
361  IF(coordinate_system%NUMBER_OF_DIMENSIONS==3) THEN
362  coordinate_convert_from_rc_sp(1)=sqrt(z(1)**2+z(2)**2+z(3)**2)
363  IF(abs(z(1))>=zero_tolerance_sp.OR.abs(z(2))>zero_tolerance_sp) THEN
364  coordinate_convert_from_rc_sp(2)=atan2(z(2),z(1))
365  ELSE
366  coordinate_convert_from_rc_sp(2)=0.0_sp
367  ENDIF
368  a1=sqrt(z(1)**2+z(2)**2)
369  IF(abs(z(3))>=zero_tolerance_sp.OR.abs(a1)>zero_tolerance_sp) THEN
370  coordinate_convert_from_rc_sp(3)=atan2(z(3),a1)
371  ELSE
372  coordinate_convert_from_rc_sp(3)=0.0_sp
373  ENDIF
374  ELSE
375  CALL flagerror("Invalid number of coordinates.",err,error,*999)
376  ENDIF
377  CASE(coordinate_prolate_spheroidal_type)
378  focus=REAL(coordinate_system%focus,sp)
379  IF(coordinate_system%NUMBER_OF_DIMENSIONS==3) THEN
380  a1=z(1)**2+z(2)**2+z(3)**2-focus**2
381  a2=sqrt(a1**2+4.0_sp*(focus**2)*(z(2)**2+z(3)**2))
382  a3=2.0_sp*focus**2
383  a4=max((a2+a1)/a3,0.0_sp)
384  a5=max((a2-a1)/a3,0.0_sp)
385  a6=sqrt(a4)
386  a7=min(sqrt(a5),1.0_sp)
387  IF(abs(a7)<=1.0_sp) THEN
388  a8=asin(a7)
389  ELSE
390  a8=0.0_sp
391  CALL flag_warning("Put A8=0 since ABS(A8)>1.",err,error,*999)
392  ENDIF
393  IF((abs(z(3))<zero_tolerance_sp).OR.(abs(a6)<zero_tolerance_sp).OR.(abs(a7)<zero_tolerance_sp)) THEN
394  a9=0.0_sp
395  ELSE
396  IF(abs(a6*a7)>zero_tolerance_sp) THEN
397  a9=z(3)/(focus*a6*a7)
398  ELSE
399  a9=0.0_sp
400  CALL flag_warning("Put A9=0 since A6*A7=0.",err,error,*999)
401  ENDIF
402  IF(a9>=1.0_sp) THEN
403  a9=REAL(pi,sp)/2.0_SP
404  ELSE IF(a9<=-1.0_sp) THEN
405  a9=-REAL(pi,sp)/2.0_SP
406  ELSE
407  a9=asin(a9)
408  ENDIF
409  ENDIF
410  coordinate_convert_from_rc_sp(1)=log(a6+sqrt(a4+1.0_sp))
411  IF(abs(z(1))>=zero_tolerance_sp) THEN
412  coordinate_convert_from_rc_sp(2)=a8
413  ELSE
414  coordinate_convert_from_rc_sp(2)=REAL(pi,sp)-A8
415  ENDIF
416  IF(abs(z(2))>=zero_tolerance_sp) THEN
417  coordinate_convert_from_rc_sp(3)=mod(a9+2.0_sp*REAL(PI,SP),2.0_SP*&
418  & REAL(PI,SP))
419  ELSE
420  coordinate_convert_from_rc_sp(3)=REAL(pi,sp)-A9
421  ENDIF
422  ELSE
423  CALL flagerror("Invalid number of coordinates.",err,error,*999)
424  ENDIF
425  CASE(coordinate_oblate_spheroidal_type)
426  CALL flagerror("Not implemented.",err,error,*999)
427  CASE DEFAULT
428  CALL flagerror("Invalid coordinate type.",err,error,*999)
429  END SELECT
430 
431  exits("COORDINATE_CONVERT_FROM_RC_SP")
432  RETURN
433 999 errorsexits("COORDINATE_CONVERT_FROM_RC_SP",err,error)
434  RETURN
435  END FUNCTION coordinate_convert_from_rc_sp
436 
437  !
438  !================================================================================================================================
439  !
440 
444  FUNCTION coordinate_convert_to_rc_dp(COORDINATE_SYSTEM,X,ERR,ERROR)
445 
446  !Argument variables
447  TYPE(coordinate_system_type), INTENT(IN) :: COORDINATE_SYSTEM
448  REAL(DP), INTENT(IN) :: X(:)
449  INTEGER(INTG), INTENT(OUT) :: ERR
450  TYPE(varying_string), INTENT(OUT) :: ERROR
451  !Function variable
452  REAL(DP) :: COORDINATE_CONVERT_TO_RC_DP(size(x,1))
453  !Local variables
454  REAL(DP) :: FOCUS
455 
456  enters("COORDINATE_CONVERT_TO_RC_DP",err,error,*999)
457 
458  coordinate_convert_to_rc_dp=0.0_dp
459 
460  IF(SIZE(x,1)<coordinate_system%NUMBER_OF_DIMENSIONS) &
461  & CALL flagerror("Size of X is less than the number of dimensions.",err,error,*999)
462 
463  SELECT CASE(coordinate_system%TYPE)
464  CASE(coordinate_rectangular_cartesian_type)
465  coordinate_convert_to_rc_dp(1:coordinate_system%NUMBER_OF_DIMENSIONS)=x(1:coordinate_system%NUMBER_OF_DIMENSIONS)
466  CASE(coordinate_cylindrical_polar_type)
467  SELECT CASE(coordinate_system%NUMBER_OF_DIMENSIONS)
468  CASE(2)
469  coordinate_convert_to_rc_dp(1)=x(1)*cos(x(2))
470  coordinate_convert_to_rc_dp(2)=x(1)*sin(x(2))
471  CASE(3)
472  coordinate_convert_to_rc_dp(1)=x(1)*cos(x(2))
473  coordinate_convert_to_rc_dp(2)=x(1)*sin(x(2))
474  coordinate_convert_to_rc_dp(3)=x(3)
475  CASE DEFAULT
476  CALL flagerror("Invalid number of coordinates.",err,error,*999)
477  END SELECT
478  CASE(coordinate_spherical_polar_type)
479  IF(coordinate_system%NUMBER_OF_DIMENSIONS==3) THEN
480  coordinate_convert_to_rc_dp(1)=x(1)*cos(x(2))*cos(x(3))
481  coordinate_convert_to_rc_dp(2)=x(1)*sin(x(2))*cos(x(3))
482  coordinate_convert_to_rc_dp(3)=x(1)*sin(x(3))
483  ELSE
484  CALL flagerror("Invalid number of coordinates.",err,error,*999)
485  ENDIF
486  CASE(coordinate_prolate_spheroidal_type)
487  IF(coordinate_system%NUMBER_OF_DIMENSIONS==3) THEN
488  focus=coordinate_system%FOCUS
489  coordinate_convert_to_rc_dp(1)=focus*cosh(x(1))*cos(x(2))
490  coordinate_convert_to_rc_dp(2)=focus*sinh(x(1))*sin(x(2))*cos(x(3))
491  coordinate_convert_to_rc_dp(3)=focus*sinh(x(1))*sin(x(2))*sin(x(3))
492  ELSE
493  CALL flagerror("Invalid number of coordinates.",err,error,*999)
494  ENDIF
495  CASE(coordinate_oblate_spheroidal_type)
496  IF(coordinate_system%NUMBER_OF_DIMENSIONS==3) THEN
497  focus=coordinate_system%FOCUS
498  coordinate_convert_to_rc_dp(1)=focus*cosh(x(1))*cos(x(2))*cos(x(3))
499  coordinate_convert_to_rc_dp(2)=focus*sinh(x(1))*sin(x(2))
500  coordinate_convert_to_rc_dp(3)=focus*cosh(x(1))*cos(x(2))*sin(x(3))
501  ELSE
502  CALL flagerror("Invalid number of coordinates.",err,error,*999)
503  ENDIF
504  CASE DEFAULT
505  CALL flagerror("Invalid coordinate type.",err,error,*999)
506  END SELECT
507 
508  exits("COORDINATE_CONVERT_TO_RC_DP")
509  RETURN
510 999 errorsexits("COORDINATE_CONVERT_TO_RC_DP",err,error)
511  RETURN
512  END FUNCTION coordinate_convert_to_rc_dp
513 
514  !
515  !================================================================================================================================
516  !
517 
521  FUNCTION coordinate_convert_to_rc_sp(COORDINATE_SYSTEM,X,ERR,ERROR)
522 
523  !Argument variables
524  TYPE(coordinate_system_type), INTENT(IN) :: COORDINATE_SYSTEM
525  REAL(SP), INTENT(IN) :: X(:)
526  INTEGER(INTG), INTENT(OUT) :: ERR
527  TYPE(varying_string), INTENT(OUT) :: ERROR
528  !Function variable
529  REAL(SP) :: COORDINATE_CONVERT_TO_RC_SP(size(x,1))
530  !Local variables
531  REAL(SP) :: FOCUS
532 
533  enters("COORDINATE_CONVERT_TO_RC_SP",err,error,*999)
534 
535  coordinate_convert_to_rc_sp=0.0_sp
536 
537  IF(SIZE(x,1)<coordinate_system%NUMBER_OF_DIMENSIONS) &
538  & CALL flagerror("Size of X is less than the number of dimensions.",err,error,*999)
539 
540  SELECT CASE(coordinate_system%TYPE)
541  CASE(coordinate_rectangular_cartesian_type)
542  coordinate_convert_to_rc_sp(1:coordinate_system%NUMBER_OF_DIMENSIONS)=x(1:coordinate_system%NUMBER_OF_DIMENSIONS)
543  CASE(coordinate_cylindrical_polar_type)
544  SELECT CASE(coordinate_system%NUMBER_OF_DIMENSIONS)
545  CASE(2)
546  coordinate_convert_to_rc_sp(1)=x(1)*cos(x(2))
547  coordinate_convert_to_rc_sp(2)=x(1)*sin(x(2))
548  CASE(3)
549  coordinate_convert_to_rc_sp(1)=x(1)*cos(x(2))
550  coordinate_convert_to_rc_sp(2)=x(1)*sin(x(2))
551  coordinate_convert_to_rc_sp(3)=x(3)
552  CASE DEFAULT
553  CALL flagerror("Invalid number of coordinates.",err,error,*999)
554  END SELECT
555  CASE(coordinate_spherical_polar_type)
556  IF(coordinate_system%NUMBER_OF_DIMENSIONS==3) THEN
557  coordinate_convert_to_rc_sp(1)=x(1)*cos(x(2))*cos(x(3))
558  coordinate_convert_to_rc_sp(2)=x(1)*sin(x(2))*cos(x(3))
559  coordinate_convert_to_rc_sp(3)=x(1)*sin(x(3))
560  ELSE
561  CALL flagerror("Invalid number of coordinates.",err,error,*999)
562  ENDIF
563  CASE(coordinate_prolate_spheroidal_type)
564  IF(coordinate_system%NUMBER_OF_DIMENSIONS==3) THEN
565  focus=REAL(coordinate_system%focus,sp)
566  coordinate_convert_to_rc_sp(1)=focus*cosh(x(1))*cos(x(2))
567  coordinate_convert_to_rc_sp(2)=focus*sinh(x(1))*sin(x(2))*cos(x(3))
568  coordinate_convert_to_rc_sp(3)=focus*sinh(x(1))*sin(x(2))*sin(x(3))
569  ELSE
570  CALL flagerror("Invalid number of coordinates.",err,error,*999)
571  ENDIF
572  CASE(coordinate_oblate_spheroidal_type)
573  IF(coordinate_system%NUMBER_OF_DIMENSIONS==3) THEN
574  focus=REAL(coordinate_system%focus,sp)
575  coordinate_convert_to_rc_sp(1)=focus*cosh(x(1))*cos(x(2))*cos(x(3))
576  coordinate_convert_to_rc_sp(2)=focus*sinh(x(1))*sin(x(2))
577  coordinate_convert_to_rc_sp(3)=focus*cosh(x(1))*cos(x(2))*sin(x(3))
578  ELSE
579  CALL flagerror("Invalid number of coordinates.",err,error,*999)
580  ENDIF
581  CASE DEFAULT
582  CALL flagerror("Invalid coordinate type.",err,error,*999)
583  END SELECT
584 
585  exits("COORDINATE_CONVERT_TO_RC_SP")
586  RETURN
587 999 errorsexits("COORDINATE_CONVERT_TO_RC_SP",err,error)
588  RETURN
589  END FUNCTION coordinate_convert_to_rc_sp
590 
591  !
592  !================================================================================================================================
593  !
594 
597  FUNCTION coordinate_delta_calculate_dp(COORDINATE_SYSTEM,X,Y,ERR,ERROR)
598 
599  !Argument variables
600  TYPE(coordinate_system_type), INTENT(IN) :: COORDINATE_SYSTEM
601  REAL(DP), INTENT(IN) :: X(:)
602  REAL(DP), INTENT(IN) :: Y(:)
603  INTEGER(INTG), INTENT(OUT) :: ERR
604  TYPE(varying_string), INTENT(OUT) :: ERROR
605  !Function variable
606  REAL(DP) :: COORDINATE_DELTA_CALCULATE_DP(size(x,1))
607  !Local variables
608 
609  enters("COORDINATE_DELTA_CALCULATE_DP",err,error,*999)
610 
611  coordinate_delta_calculate_dp=0.0_dp
612 
613  IF(SIZE(x,1)<coordinate_system%NUMBER_OF_DIMENSIONS) &
614  & CALL flagerror("Size of X is less than the number of dimensions.",err,error,*999)
615 
616  IF(SIZE(x,1)/=SIZE(y,1)) &
617  & CALL flagerror("Size of X is different to the size of Y.",err,error,*999)
618 
619  coordinate_delta_calculate_dp(1:coordinate_system%NUMBER_OF_DIMENSIONS)=y(1:coordinate_system%NUMBER_OF_DIMENSIONS)- &
620  & x(1:coordinate_system%NUMBER_OF_DIMENSIONS)
621  SELECT CASE(coordinate_system%TYPE)
622  CASE(coordinate_rectangular_cartesian_type)
623  !Do nothing
624  CASE(coordinate_cylindrical_polar_type)
625  CALL flagerror("Not implemented.",err,error,*999)
626  CASE(coordinate_spherical_polar_type)
627  CALL flagerror("Not implemented.",err,error,*999)
628  CASE(coordinate_prolate_spheroidal_type)
629  CALL flagerror("Not implemented.",err,error,*999)
630  CASE(coordinate_oblate_spheroidal_type)
631  CALL flagerror("Not implemented.",err,error,*999)
632  CASE DEFAULT
633  CALL flagerror("Invalid coordinate type.",err,error,*999)
634  END SELECT
635 
636  exits("COORDINATE_DELTA_CALCULATE_DP")
637  RETURN
638 999 errorsexits("COORDINATE_DELTA_CALCULATE_DP",err,error)
639  RETURN
640  END FUNCTION coordinate_delta_calculate_dp
641 
642  !
643  !================================================================================================================================
644  !
645 
647  SUBROUTINE coordinate_metrics_calculate(COORDINATE_SYSTEM,JACOBIAN_TYPE,METRICS,ERR,ERROR,*)
648 
649  !Argument variables
650  TYPE(coordinate_system_type), POINTER :: COORDINATE_SYSTEM
651  INTEGER(INTG), INTENT(IN) :: JACOBIAN_TYPE
652  TYPE(field_interpolated_point_metrics_type), POINTER :: METRICS
653  INTEGER(INTG), INTENT(OUT) :: ERR
654  TYPE(varying_string), INTENT(OUT) :: ERROR
655  !Local Variables
656  INTEGER(INTG) :: mi,ni,nu
657  REAL(DP) :: DET_GL,DET_DX_DXI,DX_DXI2(3),DX_DXI3(3),FF,G1,G3,LENGTH,MU,R,RC,RCRC,RR,SCALE
658  TYPE(field_interpolated_point_type), POINTER :: INTERPOLATED_POINT
659  TYPE(varying_string) :: LOCAL_ERROR
660 
661  NULLIFY(interpolated_point)
662 
663  enters("COORDINATE_METRICS_CALCULATE",err,error,*999)
664 
665  IF(ASSOCIATED(coordinate_system)) THEN
666  IF(ASSOCIATED(metrics)) THEN
667  interpolated_point=>metrics%INTERPOLATED_POINT
668  IF(ASSOCIATED(interpolated_point)) THEN
669  IF(interpolated_point%PARTIAL_DERIVATIVE_TYPE>=first_part_deriv) THEN
670 
671  SELECT CASE(metrics%NUMBER_OF_XI_DIMENSIONS)
672  CASE(1)
673  !Calculate the derivatives of X with respect to XI
674  nu=partial_derivative_first_derivative_map(1)
675  metrics%DX_DXI(1:metrics%NUMBER_OF_X_DIMENSIONS,1)=interpolated_point%VALUES(1:metrics%NUMBER_OF_X_DIMENSIONS,nu)
676  !Initialise the covariant metric tensor to the identity matrix
677  metrics%GL(1,1)=1.0_dp
678  CASE(2)
679  !Calculate the derivatives of X with respect to XI
680  nu=partial_derivative_first_derivative_map(1)
681  metrics%DX_DXI(1:metrics%NUMBER_OF_X_DIMENSIONS,1)=interpolated_point%VALUES(1:metrics%NUMBER_OF_X_DIMENSIONS,nu)
682  nu=partial_derivative_first_derivative_map(2)
683  metrics%DX_DXI(1:metrics%NUMBER_OF_X_DIMENSIONS,2)=interpolated_point%VALUES(1:metrics%NUMBER_OF_X_DIMENSIONS,nu)
684  !Initialise the covariant metric tensor to the identity matrix
685  metrics%GL(1,1)=1.0_dp
686  metrics%GL(1,2)=0.0_dp
687  metrics%GL(2,1)=0.0_dp
688  metrics%GL(2,2)=1.0_dp
689  CASE(3)
690  !Calculate the derivatives of X with respect to XI
691  nu=partial_derivative_first_derivative_map(1)
692  metrics%DX_DXI(1:metrics%NUMBER_OF_X_DIMENSIONS,1)=interpolated_point%VALUES(1:metrics%NUMBER_OF_X_DIMENSIONS,nu)
693  nu=partial_derivative_first_derivative_map(2)
694  metrics%DX_DXI(1:metrics%NUMBER_OF_X_DIMENSIONS,2)=interpolated_point%VALUES(1:metrics%NUMBER_OF_X_DIMENSIONS,nu)
695  nu=partial_derivative_first_derivative_map(3)
696  metrics%DX_DXI(1:metrics%NUMBER_OF_X_DIMENSIONS,3)=interpolated_point%VALUES(1:metrics%NUMBER_OF_X_DIMENSIONS,nu)
697  !Initialise the covariant metric tensor to the identity matrix
698  metrics%GL(1,1)=1.0_dp
699  metrics%GL(1,2)=0.0_dp
700  metrics%GL(1,3)=0.0_dp
701  metrics%GL(2,1)=0.0_dp
702  metrics%GL(2,2)=1.0_dp
703  metrics%GL(2,3)=0.0_dp
704  metrics%GL(3,1)=0.0_dp
705  metrics%GL(3,2)=0.0_dp
706  metrics%GL(3,3)=1.0_dp
707  CASE DEFAULT
708  CALL flagerror("Not implemented.",err,error,*999)
709  END SELECT
710 
711  !Calculate the covariant metric tensor GL(i,j)
712  SELECT CASE(coordinate_system%TYPE)
713  CASE(coordinate_rectangular_cartesian_type)
714  SELECT CASE(metrics%NUMBER_OF_X_DIMENSIONS)
715  CASE(1)
716  DO mi=1,metrics%NUMBER_OF_XI_DIMENSIONS
717  DO ni=1,metrics%NUMBER_OF_XI_DIMENSIONS
718  metrics%GL(mi,ni)=metrics%DX_DXI(1,mi)*metrics%DX_DXI(1,ni)
719  ENDDO !ni
720  ENDDO !mi
721  CASE(2)
722  DO mi=1,metrics%NUMBER_OF_XI_DIMENSIONS
723  DO ni=1,metrics%NUMBER_OF_XI_DIMENSIONS
724  metrics%GL(mi,ni)=metrics%DX_DXI(1,mi)*metrics%DX_DXI(1,ni)+ &
725  & metrics%DX_DXI(2,mi)*metrics%DX_DXI(2,ni)
726  ENDDO !ni
727  ENDDO !mi
728  CASE(3)
729  DO mi=1,metrics%NUMBER_OF_XI_DIMENSIONS
730  DO ni=1,metrics%NUMBER_OF_XI_DIMENSIONS
731  metrics%GL(mi,ni)=metrics%DX_DXI(1,mi)*metrics%DX_DXI(1,ni)+ &
732  & metrics%DX_DXI(2,mi)*metrics%DX_DXI(2,ni)+ &
733  & metrics%DX_DXI(3,mi)*metrics%DX_DXI(3,ni)
734  ENDDO !ni
735  ENDDO !mi
736  CASE DEFAULT
737  local_error=trim(number_to_vstring(metrics%NUMBER_OF_X_DIMENSIONS,"*",err,error))// &
738  & " is an invalid number of dimensions for a rectangular cartesian coordinate system."
739  CALL flagerror(local_error,err,error,*999)
740  END SELECT
741  CASE(coordinate_cylindrical_polar_type)
742  r=interpolated_point%VALUES(1,1)
743  rr=r*r
744  IF(metrics%NUMBER_OF_X_DIMENSIONS==2) THEN
745  DO mi=1,metrics%NUMBER_OF_XI_DIMENSIONS
746  DO ni=1,metrics%NUMBER_OF_XI_DIMENSIONS
747  metrics%GL(mi,ni)=metrics%DX_DXI(1,mi)*metrics%DX_DXI(1,ni)+rr*metrics%DX_DXI(2,mi)*metrics%DX_DXI(2,ni)
748  ENDDO !ni
749  ENDDO !mi
750  ELSE IF(metrics%NUMBER_OF_X_DIMENSIONS==3) THEN
751  DO mi=1,metrics%NUMBER_OF_XI_DIMENSIONS
752  DO ni=1,metrics%NUMBER_OF_XI_DIMENSIONS
753  metrics%GL(mi,ni)=metrics%DX_DXI(1,mi)*metrics%DX_DXI(1,ni)+rr*metrics%DX_DXI(2,mi)*metrics%DX_DXI(2,ni)+ &
754  & metrics%DX_DXI(3,mi)*metrics%DX_DXI(3,ni)
755  ENDDO !ni
756  ENDDO !mi
757  ELSE
758  local_error=trim(number_to_vstring(metrics%NUMBER_OF_X_DIMENSIONS,"*",err,error))// &
759  & " is an invalid number of dimensions for a cylindrical polar coordinate system."
760  CALL flagerror(local_error,err,error,*999)
761  ENDIF
762  CASE(coordinate_spherical_polar_type)
763  r=interpolated_point%VALUES(1,1)
764  rr=r*r
765  rc=r*cos(interpolated_point%VALUES(3,1))
766  rcrc=rc*rc
767  DO mi=1,metrics%NUMBER_OF_XI_DIMENSIONS
768  DO ni=1,metrics%NUMBER_OF_XI_DIMENSIONS
769  metrics%GL(mi,ni)=metrics%DX_DXI(1,mi)*metrics%DX_DXI(1,ni)+rcrc*metrics%DX_DXI(2,mi)*metrics%DX_DXI(2,ni)+ &
770  & rr*metrics%DX_DXI(3,mi)*metrics%DX_DXI(3,ni)
771  ENDDO !ni
772  ENDDO !mi
773  CASE(coordinate_prolate_spheroidal_type)
774  IF(abs(interpolated_point%VALUES(2,1))<zero_tolerance) THEN
775  CALL flag_warning("Mu is zero.",err,error,*999)
776  ELSE
777  ff=coordinate_system%FOCUS*coordinate_system%FOCUS
778  r=interpolated_point%VALUES(1,1)
779  mu=interpolated_point%VALUES(2,1)
780  g1=ff*(sinh(r)*sinh(r)+sin(mu)*sin(mu))
781  g3=ff*sinh(r)*sinh(r)*sin(mu)*sin(mu)
782  IF(metrics%NUMBER_OF_X_DIMENSIONS==2) THEN
783  DO mi=1,metrics%NUMBER_OF_XI_DIMENSIONS
784  DO ni=1,metrics%NUMBER_OF_XI_DIMENSIONS
785  metrics%GL(mi,ni)=g1*(metrics%DX_DXI(1,mi)*metrics%DX_DXI(1,ni)+metrics%DX_DXI(2,mi)*metrics%DX_DXI(2,ni))
786  ENDDO !ni
787  ENDDO !mi
788  ELSE IF(metrics%NUMBER_OF_X_DIMENSIONS==3) THEN
789  DO mi=1,metrics%NUMBER_OF_XI_DIMENSIONS
790  DO ni=1,metrics%NUMBER_OF_XI_DIMENSIONS
791  metrics%GL(mi,ni)=g1*(metrics%DX_DXI(1,mi)*metrics%DX_DXI(1,ni)+metrics%DX_DXI(2,mi)*metrics%DX_DXI(2,ni))+ &
792  & g3*metrics%DX_DXI(3,mi)*metrics%DX_DXI(3,ni)
793  ENDDO !ni
794  ENDDO !mi
795  ELSE
796  local_error=trim(number_to_vstring(metrics%NUMBER_OF_X_DIMENSIONS,"*",err,error))// &
797  & " is an invalid number of dimensions for a prolate spheroidal coordinate system."
798  CALL flagerror(local_error,err,error,*999)
799  ENDIF
800  ENDIF
801  CASE(coordinate_oblate_spheroidal_type)
802  CALL flagerror("Not implemented.",err,error,*999)
803  CASE DEFAULT
804  local_error="The coordinate system type of "//trim(number_to_vstring(coordinate_system%TYPE,"*",err,error))// &
805  & " is invalid."
806  CALL flagerror(local_error,err,error,*999)
807  END SELECT
808 
809  !Calcualte the contravariant metric tensor
810  CALL invert(metrics%GL(1:metrics%NUMBER_OF_XI_DIMENSIONS,1:metrics%NUMBER_OF_XI_DIMENSIONS), &
811  & metrics%GU(1:metrics%NUMBER_OF_XI_DIMENSIONS,1:metrics%NUMBER_OF_XI_DIMENSIONS),det_gl, &
812  & err,error,*999)
813 
814  !Calculate the Jacobian
815  SELECT CASE(jacobian_type)
816  CASE(coordinate_jacobian_no_type)
817  metrics%JACOBIAN=0.0
818  metrics%JACOBIAN_TYPE=coordinate_jacobian_no_type
819  CASE(coordinate_jacobian_line_type)
820  metrics%JACOBIAN=sqrt(abs(metrics%GL(1,1)))
821  metrics%JACOBIAN_TYPE=coordinate_jacobian_line_type
822  CASE(coordinate_jacobian_area_type)
823  IF(metrics%NUMBER_OF_XI_DIMENSIONS==3) THEN
824  metrics%JACOBIAN=sqrt(abs(det_gl*metrics%GU(3,3)))
825  ELSE
826  metrics%JACOBIAN=sqrt(abs(det_gl))
827  ENDIF
828  metrics%JACOBIAN_TYPE=coordinate_jacobian_area_type
829  CASE(coordinate_jacobian_volume_type)
830  metrics%JACOBIAN=sqrt(abs(det_gl))
831  metrics%JACOBIAN_TYPE=coordinate_jacobian_volume_type
832  CASE DEFAULT
833  local_error="The Jacobian type of "//trim(number_to_vstring(jacobian_type,"*",err,error))// &
834  & " is invalid."
835  CALL flagerror(local_error,err,error,*999)
836  END SELECT
837 
838  !Calculate the derivatives of Xi with respect to X - DXI_DX
839  IF(metrics%NUMBER_OF_XI_DIMENSIONS==metrics%NUMBER_OF_X_DIMENSIONS) THEN
840  CALL invert(metrics%DX_DXI,metrics%DXI_DX,det_dx_dxi,err,error,*999)
841  ELSE
842  !We have a line or a surface embedded in a higher dimensional space
843  SELECT CASE(metrics%NUMBER_OF_XI_DIMENSIONS)
844  CASE(1)
845  !Line in space
846  SELECT CASE(metrics%NUMBER_OF_X_DIMENSIONS)
847  CASE(2)
848  IF(interpolated_point%PARTIAL_DERIVATIVE_TYPE>first_part_deriv) THEN
849  !We have curvature information. Form the frenet vector frame.
850  !Calculate the normal vector from the normalised second derivative of the position vector.
851  nu=partial_derivative_second_derivative_map(1)
852  dx_dxi2=normalise(interpolated_point%VALUES(1:2,nu),err,error)
853  IF(err/=0) GOTO 999
854  ELSE
855  !No curvature information but obtain other normal frenet vector by rotating tangent vector 90 deg.
856  dx_dxi2(1)=-1.0_dp*metrics%DX_DXI(2,1)
857  dx_dxi2(2)=metrics%DX_DXI(1,1)
858  ENDIF
859  det_dx_dxi=metrics%DX_DXI(1,1)*dx_dxi2(2)-metrics%DX_DXI(2,1)*dx_dxi2(1)
860  IF(abs(det_dx_dxi)>zero_tolerance) THEN
861  metrics%DXI_DX(1,1)=dx_dxi2(2)/det_dx_dxi
862  metrics%DXI_DX(1,2)=-1.0_dp*dx_dxi2(1)/det_dx_dxi
863  !Normalise to ensure that g^11=g^1.g^1
864  length=l2norm(metrics%DXI_DX(1,1:2))
865  scale=sqrt(abs(metrics%GU(1,1)))/length
866  metrics%DXI_DX(1,1:2)=scale*metrics%DXI_DX(1,1:2)
867  ELSE
868  CALL flag_warning("Zero determinant. Unable to obtain dxi/dx.",err,error,*999)
869  metrics%DXI_DX=0.0_dp
870  ENDIF
871  CASE(3)
872  IF(interpolated_point%PARTIAL_DERIVATIVE_TYPE>first_part_deriv) THEN
873  !We have curvature information. Form the frenet vector frame.
874  !Calculate the normal vector from the normalised second derivative of the position vector.
875  nu=partial_derivative_second_derivative_map(1)
876  dx_dxi2=normalise(interpolated_point%VALUES(1:3,nu),err,error)
877  IF(err/=0) GOTO 999
878  !Calculate the bi-normal vector from the normalised cross product of the tangent and normal vectors
879  CALL norm_cross_product(metrics%DX_DXI(1:3,1),dx_dxi2,dx_dxi3,err,error,*999)
880  det_dx_dxi=metrics%DX_DXI(1,1)*(dx_dxi2(2)*dx_dxi3(3)-dx_dxi2(3)*dx_dxi3(2))+ &
881  & dx_dxi2(1)*(metrics%DX_DXI(3,1)*dx_dxi3(2)-dx_dxi3(3)*metrics%DX_DXI(2,1))+ &
882  & dx_dxi3(1)*(metrics%DX_DXI(2,1)*dx_dxi2(3)-metrics%DX_DXI(3,1)*dx_dxi2(2))
883  IF(abs(det_dx_dxi)>zero_tolerance) THEN
884  metrics%DXI_DX(1,1)=(dx_dxi3(3)*dx_dxi2(2)-dx_dxi2(3)*dx_dxi3(2))/det_dx_dxi
885  metrics%DXI_DX(1,2)=-1.0_dp*(dx_dxi3(3)*dx_dxi2(1)-dx_dxi2(3)*dx_dxi3(1))/det_dx_dxi
886  metrics%DXI_DX(1,3)=(dx_dxi3(2)*dx_dxi2(1)-dx_dxi2(2)*dx_dxi3(1))/det_dx_dxi
887  !Normalise to ensure that g^11=g^1.g^1
888  length=l2norm(metrics%DXI_DX(1,1:3))
889  scale=sqrt(abs(metrics%GU(1,1)))/length
890  metrics%DXI_DX(1,1:3)=scale*metrics%DX_DXI(1,1:3)
891  ELSE
892  CALL flag_warning("Zero determinant. Unable to obtain dxi/dx.",err,error,*999)
893  metrics%DXI_DX=0.0_dp
894  ENDIF
895  ELSE
896  metrics%DXI_DX(1,1)=metrics%DX_DXI(1,1)
897  metrics%DXI_DX(1,2)=metrics%DX_DXI(2,1)
898  metrics%DXI_DX(1,3)=metrics%DX_DXI(3,1)
899  !Normalise to ensure that g^11=g^1.g^1
900  length=l2norm(metrics%DXI_DX(1,1:3))
901  scale=sqrt(abs(metrics%GU(1,1)))/length
902  metrics%DXI_DX(1,1:3)=scale*metrics%DXI_DX(1,1:3)
903  ENDIF
904  CASE DEFAULT
905  CALL flagerror("Invalid embedding of a line in space.",err,error,*999)
906  END SELECT
907  CASE(2)
908  !Surface in space
909  IF(metrics%NUMBER_OF_X_DIMENSIONS==3) THEN
910  !Surface in 3D space.
911  !Calculate the covariant vectors g^1 and g^2. These are calculated as follows:
912  !First define g_3=g_1 x g_2, and then define g^1=((g_2 x g_3)_b)/DET_GL and g^2=((g_3 x g_1)_b)/DET_GL.
913  !The _b means lowering the index with the metric tensor of the curvilinear coordinate system.
914  !This way we have a consistent set of covariant and covariant vectors, i.e. <g_M,g^N>=delta_M^N.
915  metrics%DXI_DX(1,1:3)=(metrics%GL(2,2)*metrics%DX_DXI(1:3,1)-metrics%GL(1,2)*metrics%DX_DXI(1:3,2))/det_gl
916  metrics%DXI_DX(2,1:3)=(metrics%GL(1,1)*metrics%DX_DXI(1:3,2)-metrics%GL(2,1)*metrics%DX_DXI(1:3,1))/det_gl
917  SELECT CASE(coordinate_system%TYPE)
918  CASE(coordinate_rectangular_cartesian_type)
919  !Do nothing
920  CASE(coordinate_cylindrical_polar_type)
921  r=interpolated_point%VALUES(1,1)
922  rr=r*r
923  metrics%DXI_DX(1:2,2)=metrics%DXI_DX(1:2,2)*rr
924  CASE(coordinate_spherical_polar_type)
925  r=interpolated_point%VALUES(1,1)
926  rr=r*r
927  rc=r*cos(interpolated_point%VALUES(3,1))
928  rcrc=rc*rc
929  metrics%DXI_DX(1:2,2)=metrics%DXI_DX(1:2,2)*rcrc
930  metrics%DXI_DX(1:2,3)=metrics%DXI_DX(1:2,3)*rr
931  CASE(coordinate_prolate_spheroidal_type)
932  IF(abs(interpolated_point%VALUES(2,1))<zero_tolerance) THEN
933  CALL flag_warning("Mu is zero.",err,error,*999)
934  ELSE
935  ff=coordinate_system%FOCUS*coordinate_system%FOCUS
936  r=interpolated_point%VALUES(1,1)
937  mu=interpolated_point%VALUES(2,1)
938  g1=ff*(sinh(r)*sinh(r)+sin(mu)*sin(mu))
939  g3=ff*sinh(r)*sinh(r)*sin(mu)*sin(mu)
940  metrics%DXI_DX(1:2,1)=metrics%DXI_DX(1:2,1)*g1
941  metrics%DXI_DX(1:2,2)=metrics%DXI_DX(1:2,2)*g1
942  metrics%DXI_DX(1:2,3)=metrics%DXI_DX(1:2,3)*g3
943  ENDIF
944  CASE(coordinate_oblate_spheroidal_type)
945  CALL flagerror("Not implemented.",err,error,*999)
946  CASE DEFAULT
947  local_error="The coordinate system type of "//trim(number_to_vstring(coordinate_system%TYPE,"*",err,error))// &
948  & " is invalid."
949  CALL flagerror(local_error,err,error,*999)
950  END SELECT
951  ELSE
952  CALL flagerror("Invalid embedding of a surface in space.",err,error,*999)
953  ENDIF
954  CASE DEFAULT
955  CALL flagerror("Invalid embedding in space.",err,error,*999)
956  END SELECT
957  ENDIF
958  ELSE
959  CALL flagerror("Metrics interpolated point has not been interpolated to include first derivatives.",err,error,*999)
960  ENDIF
961  ELSE
962  CALL flagerror("Metrics interpolated point is not associated.",err,error,*999)
963  ENDIF
964  ELSE
965  CALL flagerror("Metrics is not associated.",err,error,*999)
966  ENDIF
967  ELSE
968  CALL flagerror("Coordinate system is not associated.",err,error,*999)
969  ENDIF
970 
971  IF(diagnostics1) THEN
972  CALL writestring(diagnostic_output_type,"",err,error,*999)
973  CALL write_string(diagnostic_output_type,"Coordinate system metrics:",err,error,*999)
974  CALL write_string_value(diagnostic_output_type," Coordinate system type = ",trim(coordinate_system_type_string( &
975  & coordinate_system%TYPE)),err,error,*999)
976  CALL write_string_value(diagnostic_output_type," Number of X dimensions = ",metrics%NUMBER_OF_X_DIMENSIONS,err,error,*999)
977  CALL write_string_value(diagnostic_output_type," Number of Xi dimensions = ",metrics%NUMBER_OF_XI_DIMENSIONS,err,error,*999)
978  CALL write_string(diagnostic_output_type," Location of metrics:",err,error,*999)
979  CALL write_string_vector(diagnostic_output_type,1,1,metrics%NUMBER_OF_X_DIMENSIONS,3,3,interpolated_point%VALUES(:,1), &
980  & '(" X :",3(X,E13.6))','(17X,3(X,E13.6))',err,error,*999)
981  CALL write_string(diagnostic_output_type," Derivative of X wrt Xi:",err,error,*999)
982  CALL write_string_matrix(diagnostic_output_type,1,1,metrics%NUMBER_OF_X_DIMENSIONS,1,1,metrics%NUMBER_OF_XI_DIMENSIONS, &
983  & 3,3,metrics%DX_DXI,write_string_matrix_name_and_indices,'(" dX_dXi','(",I1,",:)',' :",3(X,E13.6))', &
984  & '(17X,3(X,E13.6))',err,error,*999)
985  IF(metrics%NUMBER_OF_X_DIMENSIONS/=metrics%NUMBER_OF_XI_DIMENSIONS) THEN
986  CALL write_string(diagnostic_output_type," Constructed derivative of X wrt Xi:",err,error,*999)
987  SELECT CASE(metrics%NUMBER_OF_XI_DIMENSIONS)
988  CASE(1)
989  !Line in space
990  SELECT CASE(metrics%NUMBER_OF_X_DIMENSIONS)
991  CASE(2)
992  CALL write_string_vector(diagnostic_output_type,1,1,metrics%NUMBER_OF_X_DIMENSIONS,3,3,dx_dxi2, &
993  & '(" dX_dXi(:,2) :",3(X,E13.6))','(17X,3(X,E13.6))',err,error,*999)
994  CASE(3)
995  CALL write_string_vector(diagnostic_output_type,1,1,metrics%NUMBER_OF_X_DIMENSIONS,3,3,dx_dxi2, &
996  & '(" dX_dXi(:,2) :",3(X,E13.6))','(17X,3(X,E13.6))',err,error,*999)
997  CALL write_string_vector(diagnostic_output_type,1,1,metrics%NUMBER_OF_X_DIMENSIONS,3,3,dx_dxi3, &
998  & '(" dX_dXi(:,3) :",3(X,E13.6))','(17X,3(X,E13.6))',err,error,*999)
999  CASE DEFAULT
1000  CALL flagerror("Invalid embedding of a line in space.",err,error,*999)
1001  END SELECT
1002  CASE(2)
1003  !Surface in space
1004  SELECT CASE(metrics%NUMBER_OF_X_DIMENSIONS)
1005  CASE(3)
1006  CALL write_string_vector(diagnostic_output_type,1,1,metrics%NUMBER_OF_X_DIMENSIONS,3,3,dx_dxi3, &
1007  & '(" dX_dXi(:,3) :",3(X,E13.6))','(17X,3(X,E13.6))',err,error,*999)
1008  CASE DEFAULT
1009  CALL flagerror("Invalid embedding of a surface in space.",err,error,*999)
1010  END SELECT
1011  CASE DEFAULT
1012  CALL flagerror("Invalid embedding in space.",err,error,*999)
1013  END SELECT
1014  ENDIF
1015  CALL write_string_value(diagnostic_output_type," det dX_dXi = ",det_dx_dxi,err,error,*999)
1016  CALL write_string(diagnostic_output_type," Derivative of Xi wrt X:",err,error,*999)
1017  CALL write_string_matrix(diagnostic_output_type,1,1,metrics%NUMBER_OF_XI_DIMENSIONS,1,1,metrics%NUMBER_OF_X_DIMENSIONS, &
1018  & 3,3,metrics%DXI_DX,write_string_matrix_name_and_indices,'(" dXi_dX','(",I1,",:)',' :",3(X,E13.6))', &
1019  & '(17X,3(X,E13.6))',err,error,*999)
1020  CALL write_string(diagnostic_output_type," Covariant metric tensor:",err,error,*999)
1021  CALL write_string_matrix(diagnostic_output_type,1,1,metrics%NUMBER_OF_XI_DIMENSIONS,1,1,metrics%NUMBER_OF_XI_DIMENSIONS, &
1022  & 3,3,metrics%GL,write_string_matrix_name_and_indices,'(" GL','(",I1,",:)',' :",3(X,E13.6))','(17X,3(X,E13.6))', &
1023  & err,error,*999)
1024  CALL write_string(diagnostic_output_type," Contravariant metric tensor:",err,error,*999)
1025  CALL write_string_matrix(diagnostic_output_type,1,1,metrics%NUMBER_OF_XI_DIMENSIONS,1,1,metrics%NUMBER_OF_XI_DIMENSIONS, &
1026  & 3,3,metrics%GU,write_string_matrix_name_and_indices,'(" GU','(",I1,",:)',' :",3(X,E13.6))','(17X,3(X,E13.6))', &
1027  & err,error,*999)
1028  CALL write_string_value(diagnostic_output_type," Jacobian type = ",metrics%JACOBIAN_TYPE,err,error,*999)
1029  CALL write_string_value(diagnostic_output_type," Jacobian = ",metrics%JACOBIAN,err,error,*999)
1030  ENDIF
1031 
1032  exits("COORDINATE_METRICS_CALCULATE")
1033  RETURN
1034 999 errorsexits("COORDINATE_METRICS_CALCULATE",err,error)
1035  RETURN 1
1036  END SUBROUTINE coordinate_metrics_calculate
1037 
1038  !
1039  !================================================================================================================================
1040  !
1041 
1043  SUBROUTINE coordinate_system_normal_calculate(COORDINATE_SYSTEM,REVERSE,X,N,ERR,ERROR,*)
1044 
1045  !Argument variables
1046  TYPE(coordinate_system_type), POINTER :: COORDINATE_SYSTEM
1047  LOGICAL, INTENT(IN) :: REVERSE
1048  REAL(DP), INTENT(IN) :: X(:,:)
1049  REAL(DP), INTENT(OUT) :: N(3)
1050  INTEGER(INTG), INTENT(OUT) :: ERR
1051  TYPE(varying_string), INTENT(OUT) :: ERROR
1052  !Local Variables
1053  INTEGER(INTG) :: NUMBER_OF_X_DIMENSIONS,d_s1,d_s2,d2_s1
1054  REAL(DP) :: LENGTH,R,TANGENT1(3),TANGENT2(3)
1055  TYPE(varying_string) :: LOCAL_ERROR
1056 
1057  enters("COORDINATE_SYSTEM_NORMAL_CALCULATE",err,error,*999)
1058 
1059  IF(ASSOCIATED(coordinate_system)) THEN
1060  IF(coordinate_system%COORDINATE_SYSTEM_FINISHED) THEN
1061  CALL flagerror("Coordinate system has been finished.",err,error,*999)
1062  ELSE
1063  number_of_x_dimensions=coordinate_system%NUMBER_OF_DIMENSIONS
1064  d_s1=partial_derivative_first_derivative_map(1)
1065  d_s2=partial_derivative_first_derivative_map(2)
1066  d2_s1=partial_derivative_second_derivative_map(1)
1067  SELECT CASE(coordinate_system%TYPE)
1068  CASE(coordinate_rectangular_cartesian_type)
1069  IF(number_of_x_dimensions==2) THEN
1070  tangent1(1)=x(1,d_s1)
1071  tangent1(2)=x(2,d_s1)
1072  ELSE IF(number_of_x_dimensions==3) THEN
1073  tangent1(1)=x(1,d_s1)
1074  tangent1(2)=x(2,d_s1)
1075  tangent1(3)=x(3,d_s1)
1076  tangent2(1)=x(1,d_s2)
1077  tangent2(2)=x(2,d_s2)
1078  tangent2(3)=x(3,d_s2)
1079  ELSE
1080  local_error=trim(number_to_vstring(number_of_x_dimensions,"*",err,error))// &
1081  & " is an invalid number of dimensions to calculate a normal from in a rectangular cartesian coordinate system."
1082  CALL flagerror(local_error,err,error,*999)
1083  ENDIF
1084  CASE(coordinate_cylindrical_polar_type)
1085  r=x(1,1)
1086  IF(number_of_x_dimensions==2) THEN
1087  tangent1(1)=x(1,d_s1)*cos(x(1,d_s1))-r*sin(x(1,d_s1))*x(2,d_s1)
1088  tangent1(2)=x(2,d_s1)*sin(x(1,d_s1))+r*cos(x(1,d_s1))*x(2,d_s1)
1089  ELSE IF(number_of_x_dimensions==3) THEN
1090  tangent1(1)=x(1,d_s1)*cos(x(1,d_s1))-r*sin(x(1,d_s1))*x(2,d_s1)
1091  tangent1(2)=x(2,d_s1)*sin(x(1,d_s1))+r*cos(x(1,d_s1))*x(2,d_s1)
1092  tangent1(3)=x(3,d_s1)
1093  tangent2(1)=x(1,d_s2)*cos(x(1,d_s1))-r*sin(x(1,d_s1))*x(2,d_s2)
1094  tangent2(2)=x(1,d_s2)*sin(x(1,d_s1))+r*cos(x(1,d_s1))*x(2,d_s2)
1095  tangent2(3)=x(3,d_s2)
1096  ELSE
1097  local_error=trim(number_to_vstring(number_of_x_dimensions,"*",err,error))// &
1098  & " is an invalid number of dimensions to calculate a normal from in a rectangular cartesian coordinate system."
1099  CALL flagerror(local_error,err,error,*999)
1100  ENDIF
1101  CASE(coordinate_spherical_polar_type)
1102  r=x(1,1)
1103  tangent1(1)=x(1,d_s1)*cos(x(1,d2_s1))*cos(x(1,d_s1))- &
1104  & r*sin(x(1,d2_s1))*cos(x(1,d_s1))*x(3,d_s1)- &
1105  & r*cos(x(1,d2_s1))*sin(x(1,d_s1))*x(2,d_s1)
1106  tangent1(2)=x(1,d_s1)*cos(x(1,d2_s1))*sin(x(1,d_s1))- &
1107  & r*sin(x(1,d2_s1))*sin(x(1,d_s1))*x(3,d_s1)+ &
1108  & r*cos(x(1,d2_s1))*cos(x(1,d_s1))*x(2,d_s1)
1109  tangent1(3)=x(1,d_s1)*sin(x(1,d2_s1))+r*cos(x(1,d2_s1))*x(3,d_s1)
1110  tangent2(1)=x(1,d_s2)*cos(x(1,d2_s1))*cos(x(1,d_s1))- &
1111  & r*sin(x(1,d2_s1))*cos(x(1,d_s1))*x(3,d_s2)- &
1112  & r*cos(x(1,d2_s1))*sin(x(1,d_s1))*x(2,d_s2)
1113  tangent2(2)=x(1,d_s2)*cos(x(1,d2_s1))*sin(x(1,d_s1))- &
1114  & r*sin(x(1,d2_s1))*sin(x(1,d_s1))*x(3,d_s2)+ &
1115  & r*cos(x(1,d2_s1))*cos(x(1,d_s1))*x(2,d_s2)
1116  tangent2(3)=x(1,d_s2)*sin(x(1,d2_s1))+r*cos(x(1,d2_s1))*x(3,d_s2)
1117  CASE(coordinate_prolate_spheroidal_type)
1118  CALL flagerror("Not implemented.",err,error,*999)
1119  CASE(coordinate_oblate_spheroidal_type)
1120  CALL flagerror("Not implemented.",err,error,*999)
1121  CASE DEFAULT
1122  local_error="The coordinate system type of "//trim(number_to_vstring(coordinate_system%TYPE,"*",err,error))// &
1123  & " is invalid."
1124  CALL flagerror(local_error,err,error,*999)
1125  END SELECT
1126  IF(number_of_x_dimensions==2) THEN
1127  n(1)=-tangent1(2)
1128  n(2)=tangent1(1)
1129  length=sqrt(n(1)*n(1)+n(2)*n(2))
1130  IF(abs(length)<zero_tolerance) CALL flagerror("Zero normal vector length.",err,error,*999)
1131  IF(reverse) THEN
1132  n(1)=-n(1)/length
1133  n(2)=-n(2)/length
1134  ELSE
1135  n(1)=n(1)/length
1136  n(2)=n(2)/length
1137  ENDIF
1138  ELSE
1139  n(1)=tangent1(2)*tangent2(3)-tangent1(3)*tangent2(2)
1140  n(2)=tangent1(3)*tangent2(1)-tangent1(1)*tangent2(3)
1141  n(3)=tangent1(1)*tangent2(2)-tangent1(2)*tangent2(1)
1142  length=sqrt(n(1)*n(1)+n(2)*n(2)+n(3)*n(3))
1143  IF(reverse) THEN
1144  n(1)=-n(1)/length
1145  n(2)=-n(2)/length
1146  n(3)=-n(3)/length
1147  ELSE
1148  n(1)=n(1)/length
1149  n(2)=n(2)/length
1150  n(3)=n(3)/length
1151  ENDIF
1152  ENDIF
1153  ENDIF
1154  ELSE
1155  CALL flagerror("Coordinate system is not associated.",err,error,*999)
1156  ENDIF
1157 
1158  IF(diagnostics1) THEN
1159  CALL write_string(diagnostic_output_type,"Coordinate system metrics:",err,error,*999)
1160  CALL write_string_value(diagnostic_output_type," Coordinate system type = ",trim(coordinate_system_type_string( &
1161  & coordinate_system%TYPE)),err,error,*999)
1162  CALL write_string_value(diagnostic_output_type," Number of X dimensions = ",number_of_x_dimensions,err,error,*999)
1163  CALL write_string_vector(diagnostic_output_type,1,1,number_of_x_dimensions,3,3,x(:,1),'(" X :",3(X,E13.6))', &
1164  & '(13X,3(X,E13.6))',err,error,*999)
1165  CALL write_string_vector(diagnostic_output_type,1,1,number_of_x_dimensions,3,3,tangent1,'(" Tangent 1 :",3(X,E13.6))', &
1166  & '(13X,3(X,E13.6))',err,error,*999)
1167  IF(number_of_x_dimensions==3) THEN
1168  CALL write_string_vector(diagnostic_output_type,1,1,number_of_x_dimensions,3,3,tangent2,'(" Tangent 2 :",3(X,E13.6))', &
1169  & '(13X,3(X,E13.6))',err,error,*999)
1170  ENDIF
1171  CALL write_string_vector(diagnostic_output_type,1,1,number_of_x_dimensions,3,3,n,'(" Normal :",3(X,E13.6))', &
1172  & '(13X,3(X,E13.6))',err,error,*999)
1173  ENDIF
1174 
1175  exits("COORDINATE_SYSTEM_NORMAL_CALCULATE")
1176  RETURN
1177 999 errorsexits("COORDINATE_SYSTEM_NORMAL_CALCULATE",err,error)
1178  RETURN 1
1179  END SUBROUTINE coordinate_system_normal_calculate
1180 
1181  !
1182  !================================================================================================================================
1183  !
1184 
1186  SUBROUTINE coordinate_system_dimension_get(COORDINATE_SYSTEM,NUMBER_OF_DIMENSIONS,ERR,ERROR,*)
1187 
1188  !Argument variables
1189  TYPE(coordinate_system_type), POINTER :: COORDINATE_SYSTEM
1190  INTEGER(INTG), INTENT(OUT) :: NUMBER_OF_DIMENSIONS
1191  INTEGER(INTG), INTENT(OUT) :: ERR
1192  TYPE(varying_string), INTENT(OUT) :: ERROR
1193  !Local Variables
1194 
1195  enters("COORDINATE_SYSTEM_DIMENSION_GET",err,error,*999)
1196 
1197  IF(ASSOCIATED(coordinate_system)) THEN
1198  IF(coordinate_system%COORDINATE_SYSTEM_FINISHED) THEN
1199  number_of_dimensions=coordinate_system%NUMBER_OF_DIMENSIONS
1200  ELSE
1201  CALL flagerror("Coordinate system has not been finished.",err,error,*999)
1202  ENDIF
1203  ELSE
1204  CALL flagerror("Coordinate system is not associated.",err,error,*999)
1205  ENDIF
1206 
1207  exits("COORDINATE_SYSTEM_DIMENSION_GET")
1208  RETURN
1209 999 errorsexits("COORDINATE_SYSTEM_DIMENSION_GET",err,error)
1210  RETURN 1
1211 
1212  END SUBROUTINE coordinate_system_dimension_get
1213 
1214  !
1215  !================================================================================================================================
1216  !
1217 
1219  SUBROUTINE coordinate_system_finalise(COORDINATE_SYSTEM,ERR,ERROR,*)
1220 
1221  !Argument variables
1222  TYPE(coordinate_system_type), POINTER :: COORDINATE_SYSTEM
1223  INTEGER(INTG), INTENT(OUT) :: ERR
1224  TYPE(varying_string), INTENT(OUT) :: ERROR
1225  !Local Variables
1226 
1227  enters("COORDINATE_SYSTEM_FINALISE",err,error,*999)
1228 
1229  IF(ASSOCIATED(coordinate_system)) THEN
1230  DEALLOCATE(coordinate_system)
1231  ENDIF
1232 
1233  exits("COORDINATE_SYSTEM_FINALISE")
1234  RETURN
1235 999 errorsexits("COORDINATE_SYSTEM_FINALISE",err,error)
1236  RETURN 1
1237 
1238  END SUBROUTINE coordinate_system_finalise
1239 
1240  !
1241  !================================================================================================================================
1242  !
1243 
1245  SUBROUTINE coordinate_system_focus_get(COORDINATE_SYSTEM,FOCUS,ERR,ERROR,*)
1246 
1247  !Argument variables
1248  TYPE(coordinate_system_type), POINTER :: COORDINATE_SYSTEM
1249  REAL(DP), INTENT(OUT) :: FOCUS
1250  INTEGER(INTG), INTENT(OUT) :: ERR
1251  TYPE(varying_string), INTENT(OUT) :: ERROR
1252  !Local Variables
1253 
1254  enters("COORDINATE_SYSTEM_FOCUS_GET",err,error,*999)
1255 
1256  IF(ASSOCIATED(coordinate_system)) THEN
1257  IF(coordinate_system%COORDINATE_SYSTEM_FINISHED) THEN
1258  SELECT CASE(coordinate_system%TYPE)
1259  CASE(coordinate_prolate_spheroidal_type,coordinate_oblate_spheroidal_type)
1260  focus=coordinate_system%FOCUS
1261  CASE DEFAULT
1262  CALL flagerror("No focus defined for this coordinate system type.",err,error,*999)
1263  END SELECT
1264  ELSE
1265  CALL flagerror("Coordinate system has not been finished.",err,error,*999)
1266  ENDIF
1267  ELSE
1268  CALL flagerror("Coordinate system is not associated.",err,error,*999)
1269  ENDIF
1270 
1271  exits("COORDINATE_SYSTEM_FOCUS_GET")
1272  RETURN
1273 999 errorsexits("COORDINATE_SYSTEM_FOCUS_GET",err,error)
1274  RETURN 1
1275 
1276  END SUBROUTINE coordinate_system_focus_get
1277 
1278  !
1279  !================================================================================================================================
1280  !
1281 
1282 
1284  SUBROUTINE coordinates_radialinterpolationtypeget(coordinateSystem,radialInterpolationType,err,error,*)
1285 
1286  !Argument variables
1287  TYPE(coordinate_system_type), POINTER :: coordinateSystem
1288  INTEGER(INTG), INTENT(OUT) :: radialInterpolationType
1289  INTEGER(INTG), INTENT(OUT) :: ERR
1290  TYPE(varying_string), INTENT(OUT) :: ERROR
1291  !Local Variables
1292 
1293  enters("Coordinates_RadialInterpolationTypeGet",err,error,*999)
1294 
1295  IF(ASSOCIATED(coordinatesystem)) THEN
1296  IF(coordinatesystem%COORDINATE_SYSTEM_FINISHED) THEN
1297  SELECT CASE(coordinatesystem%TYPE)
1298  CASE(coordinate_cylindrical_polar_type,coordinate_spherical_polar_type)
1299  radialinterpolationtype=coordinatesystem%RADIAL_INTERPOLATION_TYPE
1300  CASE DEFAULT
1301  CALL flagerror("No radial interpolation type defined for this coordinate system interpolation.",err,error,*999)
1302  END SELECT
1303  ELSE
1304  CALL flagerror("Coordinate system has not been finished.",err,error,*999)
1305  ENDIF
1306  ELSE
1307  CALL flagerror("Coordinate system is not associated.",err,error,*999)
1308  ENDIF
1309 
1310  exits("Coordinates_RadialInterpolationTypeGet")
1311  RETURN
1312 999 errorsexits("Coordinates_RadialInterpolationTypeGet",err,error)
1313  RETURN 1
1314 
1315  END SUBROUTINE coordinates_radialinterpolationtypeget
1316 
1317  !
1318  !================================================================================================================================
1319  !
1320 
1322  SUBROUTINE coordinate_system_type_get(COORDINATE_SYSTEM,SYSTEM_TYPE,ERR,ERROR,*)
1323 
1324  !Argument variables
1325  TYPE(coordinate_system_type), POINTER :: COORDINATE_SYSTEM
1326  INTEGER(INTG), INTENT(OUT) :: SYSTEM_TYPE
1327  INTEGER(INTG), INTENT(OUT) :: ERR
1328  TYPE(varying_string), INTENT(OUT) :: ERROR
1329  !Local Variables
1330 
1331  enters("COORDINATE_SYSTEM_TYPE_GET",err,error,*999)
1332 
1333  IF(ASSOCIATED(coordinate_system)) THEN
1334  IF(coordinate_system%COORDINATE_SYSTEM_FINISHED) THEN
1335  system_type=coordinate_system%TYPE
1336  ELSE
1337  CALL flagerror("Coordinate system has not been finished.",err,error,*999)
1338  ENDIF
1339  ELSE
1340  CALL flagerror("Coordinate system is not associated.",err,error,*999)
1341  ENDIF
1342 
1343  exits("COORDINATE_SYSTEM_TYPE_GET")
1344  RETURN
1345 999 errorsexits("COORDINATE_SYSTEM_TYPE_GET",err,error)
1346  RETURN 1
1347 
1348  END SUBROUTINE coordinate_system_type_get
1349 
1350  !
1351  !================================================================================================================================
1352  !
1353 
1355  SUBROUTINE coordinate_system_dimension_set(COORDINATE_SYSTEM,DIMENSION,ERR,ERROR,*)
1356 
1357  !Argument variables
1358  TYPE(coordinate_system_type), POINTER :: COORDINATE_SYSTEM
1359  INTEGER(INTG), INTENT(IN) :: DIMENSION
1360  INTEGER(INTG), INTENT(OUT) :: ERR
1361  TYPE(varying_string), INTENT(OUT) :: ERROR
1362  !Local Variables
1363 
1364  enters("COORDINATE_SYSTEM_DIMENSION_SET",err,error,*999)
1365 
1366  IF(ASSOCIATED(coordinate_system)) THEN
1367  IF(coordinate_system%COORDINATE_SYSTEM_FINISHED) THEN
1368  CALL flagerror("Coordinate system has been finished.",err,error,*999)
1369  ELSE
1370  SELECT CASE(coordinate_system%TYPE)
1371  CASE(coordinate_rectangular_cartesian_type)
1372  IF(dimension>=1.AND.dimension<=3) THEN
1373  coordinate_system%NUMBER_OF_DIMENSIONS=dimension
1374  ELSE
1375  CALL flagerror("Invalid number of dimensions.",err,error,*999)
1376  ENDIF
1377  CASE(coordinate_cylindrical_polar_type)
1378  IF(dimension>=2.AND.dimension<=3) THEN
1379  coordinate_system%NUMBER_OF_DIMENSIONS=dimension
1380  ELSE
1381  CALL flagerror("Invalid number of dimensions.",err,error,*999)
1382  ENDIF
1383  CASE(coordinate_spherical_polar_type)
1384  IF(dimension==3) THEN
1385  coordinate_system%NUMBER_OF_DIMENSIONS=dimension
1386  ELSE
1387  CALL flagerror("Invalid number of dimensions.",err,error,*999)
1388  ENDIF
1389  CASE(coordinate_prolate_spheroidal_type)
1390  IF(dimension==3) THEN
1391  coordinate_system%NUMBER_OF_DIMENSIONS=dimension
1392  ELSE
1393  CALL flagerror("Invalid number of dimensions.",err,error,*999)
1394  ENDIF
1395  CASE(coordinate_oblate_spheroidal_type)
1396  IF(dimension==3) THEN
1397  coordinate_system%NUMBER_OF_DIMENSIONS=dimension
1398  ELSE
1399  CALL flagerror("Invalid number of dimensions.",err,error,*999)
1400  ENDIF
1401  CASE DEFAULT
1402  CALL flagerror("Invalid coordinate system type.",err,error,*999)
1403  END SELECT
1404  ENDIF
1405  ELSE
1406  CALL flagerror("Coordinate system is not associated.",err,error,*999)
1407  ENDIF
1408 
1409  exits("COORDINATE_SYSTEM_DIMENSION_SET")
1410  RETURN
1411 999 errorsexits("COORDINATE_SYSTEM_DIMENSION_SET",err,error)
1412  RETURN 1
1413  END SUBROUTINE coordinate_system_dimension_set
1414 
1415  !
1416  !================================================================================================================================
1417  !
1418 
1420  SUBROUTINE coordinate_system_focus_set(COORDINATE_SYSTEM,FOCUS,ERR,ERROR,*)
1421 
1422  !Argument variables
1423  TYPE(coordinate_system_type), POINTER :: COORDINATE_SYSTEM
1424  REAL(DP), INTENT(IN) :: FOCUS
1425  INTEGER(INTG), INTENT(OUT) :: ERR
1426  TYPE(varying_string), INTENT(OUT) :: ERROR
1427  !Local Variables
1428 
1429  enters("COORDINATE_SYSTEM_FOCUS_SET",err,error,*999)
1430 
1431  IF(ASSOCIATED(coordinate_system)) THEN
1432  IF(coordinate_system%COORDINATE_SYSTEM_FINISHED) THEN
1433  CALL flagerror("Coordinate system has been finished.",err,error,*999)
1434  ELSE
1435  SELECT CASE(coordinate_system%TYPE)
1436  CASE(coordinate_prolate_spheroidal_type)
1437  IF(focus>zero_tolerance) THEN
1438  coordinate_system%FOCUS=focus
1439  ELSE
1440  CALL flagerror("Focus is less than zero.",err,error,*999)
1441  ENDIF
1442  CASE(coordinate_oblate_spheroidal_type)
1443  IF(focus>zero_tolerance) THEN
1444  coordinate_system%FOCUS=focus
1445  ELSE
1446  CALL flagerror("Focus is less than zero.",err,error,*999)
1447  ENDIF
1448  CASE DEFAULT
1449  CALL flagerror("Invalid coordinate system type.",err,error,*999)
1450  END SELECT
1451  ENDIF
1452  ELSE
1453  CALL flagerror("Coordinate system is not associated.",err,error,*999)
1454  ENDIF
1455 
1456  exits("COORDINATE_SYSTEM_FOCUS_SET")
1457  RETURN
1458 999 errorsexits("COORDINATE_SYSTEM_FOCUS_SET",err,error)
1459  RETURN 1
1460  END SUBROUTINE coordinate_system_focus_set
1461 
1462  !
1463  !================================================================================================================================
1464  !
1465 
1467  SUBROUTINE coordinates_radialinterpolationtypeset(coordinateSystem,radialInterpolationType,err,error,*)
1468 
1469  !Argument variables
1470  TYPE(coordinate_system_type), POINTER :: coordinateSystem
1471  INTEGER(INTG), INTENT(IN) :: radialInterpolationType
1472  INTEGER(INTG), INTENT(OUT) :: err
1473  TYPE(varying_string), INTENT(OUT) :: error
1474  !Local Variables
1475  TYPE(varying_string) :: localError
1476 
1477  enters("Coordinates_RadialInterpolationTypeSet",err,error,*999)
1478 
1479  IF(ASSOCIATED(coordinatesystem)) THEN
1480  IF(coordinatesystem%COORDINATE_SYSTEM_FINISHED) THEN
1481  CALL flagerror("Coordinate system has been finished.",err,error,*999)
1482  ELSE
1483  SELECT CASE(coordinatesystem%TYPE)
1484  CASE(coordinate_rectangular_cartesian_type)
1485  SELECT CASE(radialinterpolationtype)
1486  CASE(coordinate_no_radial_interpolation_type)
1487  coordinatesystem%RADIAL_INTERPOLATION_TYPE=coordinate_no_radial_interpolation_type
1488  CASE DEFAULT
1489  localerror="The radial interpolation type of "//trim(number_to_vstring(radialinterpolationtype,"*",err,error))// &
1490  & " is invalid for a rectangular cartesian coordinate system."
1491  CALL flagerror(localerror,err,error,*999)
1492  END SELECT
1493  CASE(coordinate_cylindrical_polar_type,coordinate_spherical_polar_type)
1494  SELECT CASE(radialinterpolationtype)
1495  CASE(coordinate_radial_interpolation_type)
1496  coordinatesystem%RADIAL_INTERPOLATION_TYPE=coordinate_radial_interpolation_type
1497  CASE(coordinate_radial_squared_interpolation_type)
1498  coordinatesystem%RADIAL_INTERPOLATION_TYPE=coordinate_radial_squared_interpolation_type
1499  CASE DEFAULT
1500  localerror="The radial interpolation type of "//trim(number_to_vstring(radialinterpolationtype,"*",err,error))// &
1501  & " is invalid for a cylindrical/spherical coordinate system."
1502  CALL flagerror(localerror,err,error,*999)
1503  END SELECT
1504  CASE(coordinate_prolate_spheroidal_type)
1505  SELECT CASE(radialinterpolationtype)
1506  CASE(coordinate_radial_interpolation_type)
1507  coordinatesystem%RADIAL_INTERPOLATION_TYPE=coordinate_radial_interpolation_type
1508  CASE(coordinate_radial_squared_interpolation_type)
1509  coordinatesystem%RADIAL_INTERPOLATION_TYPE=coordinate_radial_squared_interpolation_type
1510  CASE(coordinate_radial_cubed_interpolation_type)
1511  coordinatesystem%RADIAL_INTERPOLATION_TYPE=coordinate_radial_cubed_interpolation_type
1512  CASE DEFAULT
1513  localerror="The radial interpolation type of "//trim(numbertovstring(radialinterpolationtype,"*",err,error))// &
1514  & " is invalid for a prolate spheroidal coordinate system."
1515  CALL flagerror(localerror,err,error,*999)
1516  END SELECT
1517  CASE(coordinate_oblate_spheroidal_type)
1518  CALL flagerror("Not implemented.",err,error,*999)
1519  CASE DEFAULT
1520  CALL flagerror("Invalid coordinate system type.",err,error,*999)
1521  END SELECT
1522  ENDIF
1523  ELSE
1524  CALL flagerror("Coordinate system is not associated.",err,error,*999)
1525  ENDIF
1526 
1527  exits("Coordinates_RadialInterpolationTypeSet")
1528  RETURN
1529 999 errorsexits("Coordinates_RadialInterpolationTypeSet",err,error)
1530  RETURN 1
1531  END SUBROUTINE coordinates_radialinterpolationtypeset
1532 
1533  !
1534  !================================================================================================================================
1535  !
1536 
1538  SUBROUTINE coordinate_system_type_set(COORDINATE_SYSTEM,TYPE,ERR,ERROR,*)
1539 
1540  !Argument variables
1541  TYPE(coordinate_system_type), POINTER :: COORDINATE_SYSTEM
1542  INTEGER(INTG), INTENT(IN) :: TYPE
1543  INTEGER(INTG), INTENT(OUT) :: ERR
1544  TYPE(varying_string), INTENT(OUT) :: ERROR
1545  !Local Variables
1546 
1547  enters("COORDINATE_SYSTEM_TYPE_SET",err,error,*999)
1548 
1549  IF(ASSOCIATED(coordinate_system)) THEN
1550  IF(coordinate_system%COORDINATE_SYSTEM_FINISHED) THEN
1551  CALL flagerror("Coordinate system has been finished.",err,error,*999)
1552  ELSE
1553  SELECT CASE(type)
1554  CASE(coordinate_rectangular_cartesian_type)
1555  coordinate_system%TYPE=coordinate_rectangular_cartesian_type
1556  CASE(coordinate_cylindrical_polar_type)
1557  coordinate_system%TYPE=coordinate_cylindrical_polar_type
1558  CASE(coordinate_spherical_polar_type)
1559  coordinate_system%TYPE=coordinate_spherical_polar_type
1560  CASE(coordinate_prolate_spheroidal_type)
1561  coordinate_system%TYPE=coordinate_prolate_spheroidal_type
1562  CASE(coordinate_oblate_spheroidal_type)
1563  coordinate_system%TYPE=coordinate_oblate_spheroidal_type
1564  CASE DEFAULT
1565  CALL flagerror("Invalid coordinate system type.",err,error,*999)
1566  END SELECT
1567  ENDIF
1568  ELSE
1569  CALL flagerror("Coordinate system is not associated.",err,error,*999)
1570  ENDIF
1571 
1572  exits("COORDINATE_SYSTEM_TYPE_SET")
1573  RETURN
1574 999 errorsexits("COORDINATE_SYSTEM_TYPE_SET",err,error)
1575  RETURN 1
1576  END SUBROUTINE coordinate_system_type_set
1577 
1578  !
1579  !================================================================================================================================
1580  !
1581 
1583  SUBROUTINE coordinate_system_origin_get(COORDINATE_SYSTEM,ORIGIN,ERR,ERROR,*)
1584 
1585  !Argument variables
1586  TYPE(coordinate_system_type), POINTER :: COORDINATE_SYSTEM
1587  REAL(DP), INTENT(OUT) :: ORIGIN(:)
1588  INTEGER(INTG), INTENT(OUT) :: ERR
1589  TYPE(varying_string), INTENT(OUT) :: ERROR
1590  !Local Variables
1591 
1592  enters("COORDINATE_SYSTEM_ORIGIN_GET",err,error,*999)
1593 
1594  IF(ASSOCIATED(coordinate_system)) THEN
1595  IF(coordinate_system%COORDINATE_SYSTEM_FINISHED) THEN
1596  IF(SIZE(origin)>=3) THEN
1597  origin(1:3)=coordinate_system%ORIGIN
1598  ELSE
1599  CALL flagerror("The origin must have >= 3 components.",err,error,*999)
1600  ENDIF
1601  ELSE
1602  CALL flagerror("Coordinate system has not been finished.",err,error,*999)
1603  ENDIF
1604  ELSE
1605  CALL flagerror("Coordinate system is not associated.",err,error,*999)
1606  ENDIF
1607 
1608  exits("COORDINATE_SYSTEM_ORIGIN_GET")
1609  RETURN
1610 999 errorsexits("COORDINATE_SYSTEM_ORIGIN_GET",err,error)
1611  RETURN 1
1612 
1613  END SUBROUTINE coordinate_system_origin_get
1614 
1615  !
1616  !================================================================================================================================
1617  !
1618 
1620  SUBROUTINE coordinate_system_origin_set(COORDINATE_SYSTEM,ORIGIN,ERR,ERROR,*)
1621 
1622  !Argument variables
1623  TYPE(coordinate_system_type), POINTER :: COORDINATE_SYSTEM
1624  REAL(DP), INTENT(IN) :: ORIGIN(:)
1625  INTEGER(INTG), INTENT(OUT) :: ERR
1626  TYPE(varying_string), INTENT(OUT) :: ERROR
1627  !Local Variables
1628 
1629  enters("COORDINATE_SYSTEM_ORIGIN_SET",err,error,*999)
1630 
1631  IF(ASSOCIATED(coordinate_system)) THEN
1632  IF(coordinate_system%COORDINATE_SYSTEM_FINISHED) THEN
1633  CALL flagerror("Coordinate system has been finished.",err,error,*999)
1634  ELSE
1635  IF(SIZE(origin)==3) THEN
1636  coordinate_system%ORIGIN=origin
1637  ELSE
1638  CALL flagerror("The origin must have exactly 3 components.",err,error,*999)
1639  ENDIF
1640  ENDIF
1641  ELSE
1642  CALL flagerror("Coordinate system is not associated.",err,error,*999)
1643  ENDIF
1644 
1645  exits("COORDINATE_SYSTEM_ORIGIN_SET")
1646  RETURN
1647 999 errorsexits("COORDINATE_SYSTEM_ORIGIN_SET",err,error)
1648  RETURN 1
1649  END SUBROUTINE coordinate_system_origin_set
1650 
1651  !
1652  !================================================================================================================================
1653  !
1654 
1656  SUBROUTINE coordinate_system_orientation_get(COORDINATE_SYSTEM,ORIENTATION,ERR,ERROR,*)
1657 
1658  !Argument variables
1659  TYPE(coordinate_system_type), POINTER :: COORDINATE_SYSTEM
1660  REAL(DP), INTENT(OUT) :: ORIENTATION(:,:)
1661  INTEGER(INTG), INTENT(OUT) :: ERR
1662  TYPE(varying_string), INTENT(OUT) :: ERROR
1663  !Local Variables
1664 
1665  enters("COORDINATE_SYSTEM_ORIENTATION_GET",err,error,*999)
1666 
1667  IF(ASSOCIATED(coordinate_system)) THEN
1668  IF(coordinate_system%COORDINATE_SYSTEM_FINISHED) THEN
1669  IF(SIZE(orientation,1)>=3.AND.SIZE(orientation,2)>=3) THEN
1670  orientation(1:3,1:3)=coordinate_system%ORIENTATION
1671  ELSE
1672  CALL flagerror("The orientation matrix must have >= 3x3 components.",err,error,*999)
1673  ENDIF
1674  ELSE
1675  CALL flagerror("Coordinate system has not been finished.",err,error,*999)
1676  ENDIF
1677  ELSE
1678  CALL flagerror("Coordinate system is not associated.",err,error,*999)
1679  ENDIF
1680 
1681  exits("COORDINATE_SYSTEM_ORIENTATION_GET")
1682  RETURN
1683 999 errorsexits("COORDINATE_SYSTEM_ORIENTATION_GET",err,error)
1684  RETURN 1
1685  END SUBROUTINE coordinate_system_orientation_get
1686 
1687  !
1688  !================================================================================================================================
1689  !
1690 
1692  SUBROUTINE coordinate_system_orientation_set(COORDINATE_SYSTEM,ORIENTATION,ERR,ERROR,*)
1693 
1694  !Argument variables
1695  TYPE(coordinate_system_type), POINTER :: COORDINATE_SYSTEM
1696  REAL(DP), INTENT(IN) :: ORIENTATION(:,:)
1697  INTEGER(INTG), INTENT(OUT) :: ERR
1698  TYPE(varying_string), INTENT(OUT) :: ERROR
1699  !Local Variables
1700 
1701  enters("COORDINATE_SYSTEM_ORIENTATION_SET",err,error,*999)
1702 
1703  IF(ASSOCIATED(coordinate_system)) THEN
1704  IF(coordinate_system%COORDINATE_SYSTEM_FINISHED) THEN
1705  CALL flagerror("Coordinate system has been finished.",err,error,*999)
1706  ELSE
1707  IF(SIZE(orientation,1)==3.AND.SIZE(orientation,2)==3) THEN
1708 !!TODO: \todo Check orientation matrix vectors are orthogonal to each other etc.
1709  coordinate_system%ORIENTATION=orientation
1710  ELSE
1711  CALL flagerror("The orientation matrix must have exactly 3x3 components.",err,error,*999)
1712  ENDIF
1713  ENDIF
1714  ELSE
1715  CALL flagerror("Coordinate system is not associated.",err,error,*999)
1716  ENDIF
1717 
1718  exits("COORDINATE_SYSTEM_ORIENTATION_SET")
1719  RETURN
1720 999 errorsexits("COORDINATE_SYSTEM_ORIENTATION_SET",err,error)
1721  RETURN 1
1722  END SUBROUTINE coordinate_system_orientation_set
1723 
1724  !
1725  !================================================================================================================================
1726  !
1727 
1736  SUBROUTINE coordinate_system_create_start(USER_NUMBER,COORDINATE_SYSTEM,ERR,ERROR,*)
1737 
1738  !Argument variables
1739  INTEGER(INTG), INTENT(IN) :: USER_NUMBER
1740  TYPE(coordinate_system_type), POINTER :: COORDINATE_SYSTEM
1741  INTEGER(INTG), INTENT(OUT) :: ERR
1742  TYPE(varying_string), INTENT(OUT) :: ERROR
1743  !Local Variables
1744  INTEGER(INTG) :: coord_system_idx
1745  TYPE(coordinate_system_type), POINTER :: NEW_COORDINATE_SYSTEM
1746  TYPE(coordinate_system_ptr_type), POINTER :: NEW_COORDINATE_SYSTEMS(:)
1747  TYPE(varying_string) :: LOCAL_ERROR
1748 
1749  NULLIFY(new_coordinate_system)
1750  NULLIFY(new_coordinate_systems)
1751 
1752  enters("COORDINATE_SYSTEM_CREATE_START",err,error,*998)
1753 
1754  NULLIFY(new_coordinate_system)
1755  CALL coordinate_system_user_number_find(user_number,new_coordinate_system,err,error,*999)
1756  IF(ASSOCIATED(new_coordinate_system)) THEN
1757  local_error="Coordinate system number "//trim(number_to_vstring(user_number,"*",err,error))// &
1758  & " has already been created."
1759  CALL flagerror(local_error,err,error,*998)
1760  ELSE
1761  IF(ASSOCIATED(coordinate_system)) THEN
1762  CALL flagerror("Coordinate system is already associated.",err,error,*999)
1763  ELSE
1764  NULLIFY(new_coordinate_system)
1765  ALLOCATE(new_coordinate_system,stat=err)
1766  IF(err/=0) CALL flagerror("Could not allocate new coordinate system.",err,error,*999)
1767 
1768  new_coordinate_system%USER_NUMBER=user_number
1769  new_coordinate_system%COORDINATE_SYSTEM_FINISHED=.false.
1770  new_coordinate_system%TYPE=coordinate_rectangular_cartesian_type
1771  new_coordinate_system%RADIAL_INTERPOLATION_TYPE=coordinate_no_radial_interpolation_type
1772  new_coordinate_system%NUMBER_OF_DIMENSIONS=3
1773  new_coordinate_system%FOCUS=1.0_dp
1774  new_coordinate_system%ORIGIN=(/0.0_dp,0.0_dp,0.0_dp/)
1775  new_coordinate_system%ORIENTATION=reshape(&
1776  & (/1.0_dp,0.0_dp,0.0_dp, &
1777  & 0.0_dp,1.0_dp,0.0_dp, &
1778  & 0.0_dp,0.0_dp,1.0_dp/), &
1779  & (/3,3/))
1780 
1781  ALLOCATE(new_coordinate_systems(coordinate_systems%NUMBER_OF_COORDINATE_SYSTEMS+1),stat=err)
1782  IF(err/=0) CALL flagerror("Could not allocate new coordinate systems.",err,error,*999)
1783  DO coord_system_idx=1,coordinate_systems%NUMBER_OF_COORDINATE_SYSTEMS
1784  new_coordinate_systems(coord_system_idx)%PTR=>coordinate_systems%COORDINATE_SYSTEMS(coord_system_idx)%PTR
1785  ENDDO !coord_system_idx
1786  new_coordinate_systems(coordinate_systems%NUMBER_OF_COORDINATE_SYSTEMS+1)%PTR=>new_coordinate_system
1787  DEALLOCATE(coordinate_systems%COORDINATE_SYSTEMS)
1788  coordinate_systems%COORDINATE_SYSTEMS=>new_coordinate_systems
1789  coordinate_systems%NUMBER_OF_COORDINATE_SYSTEMS=coordinate_systems%NUMBER_OF_COORDINATE_SYSTEMS+1
1790 
1791  coordinate_system=>new_coordinate_system
1792  ENDIF
1793  ENDIF
1794 
1795  exits("COORDINATE_SYSTEM_CREATE_START")
1796  RETURN
1797 999 IF(ASSOCIATED(new_coordinate_system)) DEALLOCATE(new_coordinate_system)
1798  IF(ASSOCIATED(new_coordinate_systems)) DEALLOCATE(new_coordinate_systems)
1799  NULLIFY(coordinate_system)
1800 998 errorsexits("COORDINATE_SYSTEM_CREATE_START",err,error)
1801  RETURN 1
1802  END SUBROUTINE coordinate_system_create_start
1803 
1804  !
1805  !================================================================================================================================
1806  !
1807 
1809  SUBROUTINE coordinate_system_create_finish(COORDINATE_SYSTEM,ERR,ERROR,*)
1810 
1811  !Argument variables
1812  TYPE(coordinate_system_type), POINTER :: COORDINATE_SYSTEM
1813  INTEGER(INTG), INTENT(OUT) :: ERR
1814  TYPE(varying_string), INTENT(OUT) :: ERROR
1815  !Local Variables
1816  INTEGER(INTG) :: coord_system_idx
1817 
1818  enters("COORDINATE_SYSTEM_CREATE_FINISH",err,error,*999)
1819 
1820  IF(ASSOCIATED(coordinate_system)) THEN
1821  coordinate_system%COORDINATE_SYSTEM_FINISHED=.true.
1822  ELSE
1823  CALL flagerror("Coordinate system is not associated.",err,error,*999)
1824  ENDIF
1825 
1826  IF(diagnostics1) THEN
1827  CALL write_string_value(diagnostic_output_type,"Number of coordinate systems = ", &
1828  & coordinate_systems%NUMBER_OF_COORDINATE_SYSTEMS,err,error,*999)
1829  DO coord_system_idx=1,coordinate_systems%NUMBER_OF_COORDINATE_SYSTEMS
1830  CALL write_string_value(diagnostic_output_type," Coordinate system : ",coord_system_idx,err,error,*999)
1831  CALL write_string_value(diagnostic_output_type," Number = ", &
1832  & coordinate_systems%COORDINATE_SYSTEMS(coord_system_idx)%PTR%USER_NUMBER,err,error,*999)
1833  CALL write_string_value(diagnostic_output_type," Type = ", &
1834  & coordinate_system_type_string(coordinate_systems%COORDINATE_SYSTEMS(coord_system_idx)%PTR%TYPE),err,error,*999)
1835  CALL write_string_value(diagnostic_output_type," Number of dimensions = ", &
1836  & coordinate_systems%COORDINATE_SYSTEMS(coord_system_idx)%PTR%NUMBER_OF_DIMENSIONS,err,error,*999)
1837  ENDDO !coord_system_idx
1838  ENDIF
1839 
1840  exits("COORDINATE_SYSTEM_CREATE_FINISH")
1841  RETURN
1842 999 errorsexits("COORDINATE_SYSTEM_CREATE_FINISH",err,error)
1843  RETURN 1
1844  END SUBROUTINE coordinate_system_create_finish
1845 
1846  !
1847  !================================================================================================================================
1848  !
1849 
1851  SUBROUTINE coordinate_system_destroy(COORDINATE_SYSTEM,ERR,ERROR,*)
1852 
1853  !Argument variables
1854  TYPE(coordinate_system_type), POINTER :: COORDINATE_SYSTEM
1855  INTEGER(INTG), INTENT(OUT) :: ERR
1856  TYPE(varying_string), INTENT(OUT) :: ERROR
1857  !Local Variables
1858  INTEGER(INTG) :: coord_system_no,new_coord_system_no
1859  LOGICAL :: FOUND
1860  TYPE(coordinate_system_ptr_type), POINTER :: NEW_COORDINATE_SYSTEMS(:)
1861 
1862  enters("COORDINATE_SYSTEM_DESTROY",err,error,*999)
1863 
1864  IF(ASSOCIATED(coordinate_system)) THEN
1865  IF(coordinate_system%USER_NUMBER==0) THEN
1866  CALL flagerror("Cannot destroy the world coordinate system.",err,error,*999)
1867  ELSE
1868  found=.false.
1869  new_coord_system_no=0
1870  ALLOCATE(new_coordinate_systems(coordinate_systems%NUMBER_OF_COORDINATE_SYSTEMS-1),stat=err)
1871  IF(err/=0) CALL flagerror("Could not allocate new coordianate systems.",err,error,*999)
1872  DO coord_system_no=1,coordinate_systems%NUMBER_OF_COORDINATE_SYSTEMS
1873  IF(coordinate_systems%COORDINATE_SYSTEMS(coord_system_no)%PTR%USER_NUMBER==coordinate_system%USER_NUMBER) THEN
1874  found=.true.
1875  ELSE
1876  new_coord_system_no=new_coord_system_no+1
1877  new_coordinate_systems(new_coord_system_no)%PTR=>coordinate_systems%COORDINATE_SYSTEMS(coord_system_no)%PTR
1878  ENDIF
1879  ENDDO !coord_system_no
1880  IF(found) THEN
1881  CALL coordinate_system_finalise(coordinate_system,err,error,*999)
1882  DEALLOCATE(coordinate_systems%COORDINATE_SYSTEMS)
1883  coordinate_systems%COORDINATE_SYSTEMS=>new_coordinate_systems
1884  coordinate_systems%NUMBER_OF_COORDINATE_SYSTEMS=coordinate_systems%NUMBER_OF_COORDINATE_SYSTEMS-1
1885  ELSE
1886  DEALLOCATE(new_coordinate_systems)
1887  CALL flagerror("Coordinate system number to destroy does not exist.",err,error,*999)
1888  ENDIF
1889  ENDIF
1890  ELSE
1891  CALL flagerror("Coordinate system is not associated.",err,error,*999)
1892  ENDIF
1893 
1894  exits("COORDINATE_SYSTEM_DESTROY")
1895  RETURN
1896 999 errorsexits("COORDINATE_SYSTEM_DESTROY",err,error)
1897  RETURN 1
1898  END SUBROUTINE coordinate_system_destroy
1899 
1900  !
1901  !================================================================================================================================
1902  !
1903 
1905  FUNCTION dxz_dp(COORDINATE_SYSTEM,I,X,ERR,ERROR)
1906 
1907  !Argument variables
1908  TYPE(coordinate_system_type), INTENT(IN) :: COORDINATE_SYSTEM
1909  INTEGER(INTG), INTENT(IN) :: I
1910  REAL(DP), INTENT(IN) :: X(:)
1911  INTEGER(INTG), INTENT(OUT) :: ERR
1912  TYPE(varying_string), INTENT(OUT) :: ERROR
1913  !Function variable
1914  REAL(DP) :: DXZ_DP(size(x,1))
1915  !Local variables
1916  REAL(DP) :: RD,FOCUS
1917 
1918  enters("DXZ_DP",err,error,*999)
1919 
1920  dxz_dp=0.0_dp
1921 
1922  IF(SIZE(x,1)<coordinate_system%NUMBER_OF_DIMENSIONS) &
1923  & CALL flagerror("Size of X is less than the number of dimensions", &
1924  & err,error,*999)
1925 
1926  SELECT CASE(coordinate_system%TYPE)
1927  CASE(coordinate_rectangular_cartesian_type)
1928  IF(i>0.AND.i<=coordinate_system%NUMBER_OF_DIMENSIONS) THEN
1929  dxz_dp(i)=1.0_dp
1930  ELSE
1931  CALL flagerror("Invalid i value",err,error,*999)
1932  ENDIF
1933  CASE(coordinate_cylindrical_polar_type)
1934  SELECT CASE(coordinate_system%NUMBER_OF_DIMENSIONS)
1935  CASE(2)
1936  SELECT CASE(i)
1937  CASE(1)
1938  dxz_dp(1)=cos(x(2))
1939  dxz_dp(2)=-sin(x(2))/x(1)
1940  CASE(2)
1941  dxz_dp(1)=sin(x(2))
1942  dxz_dp(2)=cos(x(2))/x(1)
1943  CASE DEFAULT
1944  CALL flagerror("Invalid i value",err,error,*999)
1945  END SELECT
1946  CASE(3)
1947  SELECT CASE(i)
1948  CASE(1)
1949  dxz_dp(1)=cos(x(2))
1950  dxz_dp(2)=-sin(x(2))/x(1)
1951  dxz_dp(3)=0.0_dp
1952  CASE(2)
1953  dxz_dp(1)=sin(x(2))
1954  dxz_dp(2)=cos(x(2))/x(1)
1955  dxz_dp(3)=0.0_dp
1956  CASE(3)
1957  dxz_dp(1)=0.0_dp
1958  dxz_dp(2)=0.0_dp
1959  dxz_dp(3)=1.0_dp
1960  CASE DEFAULT
1961  CALL flagerror("Invalid i value",err,error,*999)
1962  END SELECT
1963  CASE DEFAULT
1964  CALL flagerror("Invalid number of coordinates",err,error,*999)
1965  END SELECT
1966  CASE(coordinate_spherical_polar_type)
1967  IF(coordinate_system%NUMBER_OF_DIMENSIONS==3) THEN
1968  SELECT CASE(i)
1969  CASE(1)
1970  dxz_dp(1)=cos(x(2))*cos(x(3))
1971  dxz_dp(2)=-sin(x(2))/(x(1)*cos(x(3)))
1972  dxz_dp(3)=-cos(x(2))*sin(x(3))/x(1)
1973  CASE(2)
1974  dxz_dp(1)=sin(x(2))*cos(x(3))
1975  dxz_dp(2)=cos(x(2))/(x(1)*cos(x(3)))
1976  dxz_dp(3)=-sin(x(2))*sin(x(3))/x(1)
1977  CASE(3)
1978  dxz_dp(1)=sin(x(3))
1979  dxz_dp(2)=0.0_dp
1980  dxz_dp(3)=cos(x(3))/x(1)
1981  CASE DEFAULT
1982  CALL flagerror("Invalid i value",err,error,*999)
1983  END SELECT
1984  ELSE
1985  CALL flagerror("Invalid number of coordinates",err,error,*999)
1986  ENDIF
1987  CASE(coordinate_prolate_spheroidal_type)
1988  IF(coordinate_system%NUMBER_OF_DIMENSIONS==3) THEN
1989  focus=coordinate_system%FOCUS
1990  rd=focus*(cosh(x(1))*cosh(x(1))-cos(x(2))*cos(x(2)))
1991  SELECT CASE(i)
1992  CASE(1)
1993  dxz_dp(1)=sinh(x(1))*cos(x(2))/rd
1994  dxz_dp(2)=-cosh(x(1))*sin(x(2))/rd
1995  dxz_dp(3)=0.0_dp
1996  CASE(2)
1997  dxz_dp(1)=cosh(x(1))*sin(x(2))*cos(x(3))/rd
1998  dxz_dp(2)=sinh(x(1))*cos(x(2))*cos(x(3))/rd
1999  dxz_dp(3)=-sin(x(3))/(focus*sinh(x(1))*sin(x(2)))
2000  CASE(3)
2001  dxz_dp(1)=cosh(x(1))*sin(x(2))*sin(x(3))/rd
2002  dxz_dp(2)=sinh(x(1))*cos(x(2))*sin(x(3))/rd
2003  dxz_dp(3)=cos(x(3))/(focus*sinh(x(1))*sin(x(2)))
2004  CASE DEFAULT
2005  CALL flagerror("Invalid i value",err,error,*999)
2006  END SELECT
2007  ELSE
2008  CALL flagerror("Invalid number of coordinates",err,error,*999)
2009  ENDIF
2010  CASE(coordinate_oblate_spheroidal_type)
2011  CALL flagerror("Not implemented",err,error,*999)
2012  CASE DEFAULT
2013  CALL flagerror("Invalid coordinate type",err,error,*999)
2014  END SELECT
2015 
2016  exits("DXZ_DP")
2017  RETURN
2018 999 errorsexits("DXZ_DP",err,error)
2019  RETURN
2020  END FUNCTION dxz_dp
2021 
2022  !
2023  !================================================================================================================================
2024  !
2025 
2026  !#### Generic-Function: D2ZX
2027  !### Description:
2028  !### Calculates D2Z(:)/DX(I)DX(J) at X(:), where Z(:) are rectangalar
2029  !### Cartesian and X(I) and X(J) are curvilinear coordinates defined by
2030  !### COORDINATE_SYSTEM.
2031  !### Child-functions: D2ZX_SP,D2ZX_SP
2032 
2033  !
2034  !================================================================================================================================
2035  !
2036 
2037  FUNCTION d2zx_dp(COORDINATE_SYSTEM,I,J,X,ERR,ERROR)
2038 
2039  !#### Function: D2ZX_DP
2040  !### Type: REAL(DP)(SIZE(X,1))
2041  !### Description:
2042  !### Calculates D2Z(:)/DX(I)DX(J) at X(:), where Z(:) are rectangalar
2043  !### Cartesian and X(I) and X(J) are curvilinear coordinates defined by
2044  !### COORDINATE_SYSTEM.
2045  !### Parent-function: D2ZX
2046 
2047  !Argument variables
2048  TYPE(coordinate_system_type), INTENT(IN) :: COORDINATE_SYSTEM
2049  INTEGER(INTG), INTENT(IN) :: I,J
2050  REAL(DP), INTENT(IN) :: X(:)
2051  INTEGER(INTG), INTENT(OUT) :: ERR
2052  TYPE(varying_string), INTENT(OUT) :: ERROR
2053  !Function variable
2054  REAL(DP) :: D2ZX_DP(size(x,1))
2055  !Local variables
2056  REAL(DP) :: FOCUS
2057 
2058  enters("D2ZX_DP",err,error,*999)
2059 
2060  d2zx_dp=0.0_dp
2061 
2062  IF(SIZE(x,1)<coordinate_system%NUMBER_OF_DIMENSIONS) &
2063  & CALL flagerror("Size of X is less than the number of dimensions", &
2064  & err,error,*999)
2065 
2066  SELECT CASE(coordinate_system%TYPE)
2067  CASE(coordinate_rectangular_cartesian_type)
2068  d2zx_dp(1:coordinate_system%NUMBER_OF_DIMENSIONS)=0.0_dp
2069  CASE(coordinate_cylindrical_polar_type)
2070  SELECT CASE(coordinate_system%NUMBER_OF_DIMENSIONS)
2071  CASE(2)
2072  SELECT CASE(i)
2073  CASE(1)
2074  SELECT CASE(j)
2075  CASE(1)
2076  d2zx_dp(1)=0.0_dp
2077  d2zx_dp(2)=0.0_dp
2078  CASE(2)
2079  d2zx_dp(1)=-sin(x(2))
2080  d2zx_dp(2)=cos(x(2))
2081  CASE DEFAULT
2082  CALL flagerror("Invalid j value",err,error,*999)
2083  END SELECT
2084  CASE(2)
2085  SELECT CASE(j)
2086  CASE(1)
2087  d2zx_dp(1)=-sin(x(2))
2088  d2zx_dp(2)=cos(x(2))
2089  CASE(2)
2090  d2zx_dp(1)=-x(1)*cos(x(2))
2091  d2zx_dp(2)=-x(1)*sin(x(2))
2092  CASE DEFAULT
2093  CALL flagerror("Invalid j value",err,error,*999)
2094  END SELECT
2095  CASE DEFAULT
2096  CALL flagerror("Invalid i value",err,error,*999)
2097  END SELECT
2098  CASE(3)
2099  SELECT CASE(i)
2100  CASE(1)
2101  SELECT CASE(j)
2102  CASE(1)
2103  d2zx_dp(1)=0.0_dp
2104  d2zx_dp(2)=0.0_dp
2105  d2zx_dp(3)=0.0_dp
2106  CASE(2)
2107  d2zx_dp(1)=-sin(x(2))
2108  d2zx_dp(2)=cos(x(2))
2109  d2zx_dp(3)=0.0_dp
2110  CASE(3)
2111  d2zx_dp(1)=0.0_dp
2112  d2zx_dp(2)=0.0_dp
2113  d2zx_dp(3)=0.0_dp
2114  CASE DEFAULT
2115  CALL flagerror("Invalid j value",err,error,*999)
2116  END SELECT
2117  CASE(2)
2118  SELECT CASE(j)
2119  CASE(1)
2120  d2zx_dp(1)=-sin(x(2))
2121  d2zx_dp(2)=cos(x(2))
2122  d2zx_dp(3)=0.0_dp
2123  CASE(2)
2124  d2zx_dp(1)=-x(1)*cos(x(2))
2125  d2zx_dp(2)=-x(1)*sin(x(2))
2126  d2zx_dp(3)=0.0_dp
2127  CASE(3)
2128  d2zx_dp(1)=0.0_dp
2129  d2zx_dp(2)=0.0_dp
2130  d2zx_dp(3)=0.0_dp
2131  CASE DEFAULT
2132  CALL flagerror("Invalid j value",err,error,*999)
2133  END SELECT
2134  CASE(3)
2135  SELECT CASE(j)
2136  CASE(1)
2137  d2zx_dp(1)=0.0_dp
2138  d2zx_dp(2)=0.0_dp
2139  d2zx_dp(3)=0.0_dp
2140  CASE(2)
2141  d2zx_dp(1)=0.0_dp
2142  d2zx_dp(2)=0.0_dp
2143  d2zx_dp(3)=0.0_dp
2144  CASE(3)
2145  d2zx_dp(1)=0.0_dp
2146  d2zx_dp(2)=0.0_dp
2147  d2zx_dp(3)=0.0_dp
2148  CASE DEFAULT
2149  CALL flagerror("Invalid j value",err,error,*999)
2150  END SELECT
2151  CASE DEFAULT
2152  CALL flagerror("Invalid i value",err,error,*999)
2153  END SELECT
2154  CASE DEFAULT
2155  CALL flagerror("Invalid number of coordinates",err,error,*999)
2156  END SELECT
2157  CASE(coordinate_spherical_polar_type)
2158  IF(coordinate_system%NUMBER_OF_DIMENSIONS==3) THEN
2159  SELECT CASE(i)
2160  CASE(1)
2161  SELECT CASE(j)
2162  CASE(1)
2163  d2zx_dp(1)=0.0_dp
2164  d2zx_dp(2)=0.0_dp
2165  d2zx_dp(3)=0.0_dp
2166  CASE(2)
2167  d2zx_dp(1)=-sin(x(2))*cos(x(3))
2168  d2zx_dp(2)=cos(x(2))*cos(x(3))
2169  d2zx_dp(3)=0.0_dp
2170  CASE(3)
2171  d2zx_dp(1)=-cos(x(2))*sin(x(3))
2172  d2zx_dp(2)=-sin(x(2))*sin(x(3))
2173  d2zx_dp(3)=cos(x(3))
2174  CASE DEFAULT
2175  CALL flagerror("Invalid j value",err,error,*999)
2176  END SELECT
2177  CASE(2)
2178  SELECT CASE(j)
2179  CASE(1)
2180  d2zx_dp(1)=-sin(x(2))*cos(x(3))
2181  d2zx_dp(2)=cos(x(2))*cos(x(3))
2182  d2zx_dp(3)=0.0_dp
2183  CASE(2)
2184  d2zx_dp(1)=-x(1)*cos(x(2))*cos(x(3))
2185  d2zx_dp(2)=-x(1)*sin(x(2))*cos(x(3))
2186  d2zx_dp(3)=0.0_dp
2187  CASE(3)
2188  d2zx_dp(1)=x(1)*sin(x(2))*sin(x(3))
2189  d2zx_dp(2)=-x(1)*cos(x(2))*sin(x(3))
2190  d2zx_dp(3)=0.0_dp
2191  CASE DEFAULT
2192  CALL flagerror("Invalid j value",err,error,*999)
2193  END SELECT
2194  CASE(3)
2195  SELECT CASE(j)
2196  CASE(1)
2197  d2zx_dp(1)=-cos(x(2))*sin(x(3))
2198  d2zx_dp(2)=-sin(x(2))*sin(x(3))
2199  d2zx_dp(3)=cos(x(3))
2200  CASE(2)
2201  d2zx_dp(1)=x(1)*sin(x(2))*sin(x(3))
2202  d2zx_dp(2)=-x(1)*cos(x(2))*sin(x(3))
2203  d2zx_dp(3)=0.0_dp
2204  CASE(3)
2205  d2zx_dp(1)=-x(1)*cos(x(2))*cos(x(3))
2206  d2zx_dp(2)=-x(1)*sin(x(2))*cos(x(3))
2207  d2zx_dp(3)=-x(1)*sin(x(3))
2208  CASE DEFAULT
2209  CALL flagerror("Invalid j value",err,error,*999)
2210  END SELECT
2211  CASE DEFAULT
2212  CALL flagerror("Invalid i value",err,error,*999)
2213  END SELECT
2214  ELSE
2215  CALL flagerror("Invalid number of coordinates",err,error,*999)
2216  ENDIF
2217  CASE(coordinate_prolate_spheroidal_type)
2218  IF(coordinate_system%NUMBER_OF_DIMENSIONS==3) THEN
2219  focus=coordinate_system%FOCUS
2220  SELECT CASE(i)
2221  CASE(1)
2222  SELECT CASE(j)
2223  CASE(1)
2224  d2zx_dp(1)=focus*cosh(x(1))*cos(x(2))
2225  d2zx_dp(2)=focus*sinh(x(1))*sin(x(2))*cos(x(3))
2226  d2zx_dp(3)=focus*sinh(x(1))*sin(x(2))*sin(x(3))
2227  CASE(2)
2228  d2zx_dp(1)=-focus*sinh(x(1))*sin(x(2))
2229  d2zx_dp(2)=focus*cosh(x(1))*cos(x(2))*cos(x(3))
2230  d2zx_dp(3)=focus*cosh(x(1))*cos(x(2))*sin(x(3))
2231  CASE(3)
2232  d2zx_dp(1)=0.0_dp
2233  d2zx_dp(2)=-focus*cosh(x(1))*sin(x(2))*sin(x(3))
2234  d2zx_dp(3)=focus*cosh(x(1))*sin(x(2))*cos(x(3))
2235  CASE DEFAULT
2236  CALL flagerror("Invalid j value",err,error,*999)
2237  END SELECT
2238  CASE(2)
2239  SELECT CASE(j)
2240  CASE(1)
2241  d2zx_dp(1)=-focus*sinh(x(1))*sin(x(2))
2242  d2zx_dp(2)=focus*cosh(x(1))*cos(x(2))*cos(x(3))
2243  d2zx_dp(3)=focus*cosh(x(1))*cos(x(2))*sin(x(3))
2244  CASE(2)
2245  d2zx_dp(1)=-focus*cosh(x(1))*cos(x(2))
2246  d2zx_dp(2)=-focus*sinh(x(1))*sin(x(2))*cos(x(3))
2247  d2zx_dp(3)=-focus*sinh(x(1))*sin(x(2))*sin(x(3))
2248  CASE(3)
2249  d2zx_dp(1)=0.0_dp
2250  d2zx_dp(2)=-focus*sinh(x(1))*cos(x(2))*sin(x(3))
2251  d2zx_dp(3)=focus*sinh(x(1))*cos(x(2))*cos(x(3))
2252  CASE DEFAULT
2253  CALL flagerror("Invalid j value",err,error,*999)
2254  END SELECT
2255  CASE(3)
2256  SELECT CASE(j)
2257  CASE(1)
2258  d2zx_dp(1)=0.0_dp
2259  d2zx_dp(2)=-focus*cosh(x(1))*sin(x(2))*sin(x(3))
2260  d2zx_dp(3)=focus*cosh(x(1))*sin(x(2))*cos(x(3))
2261  CASE(2)
2262  d2zx_dp(1)=0.0_dp
2263  d2zx_dp(2)=-focus*sinh(x(1))*cos(x(2))*sin(x(3))
2264  d2zx_dp(3)=focus*sinh(x(1))*cos(x(2))*cos(x(3))
2265  CASE(3)
2266  d2zx_dp(1)=0.0_dp
2267  d2zx_dp(2)=-focus*sinh(x(1))*sin(x(2))*cos(x(3))
2268  d2zx_dp(3)=-focus*sinh(x(1))*sin(x(2))*sin(x(3))
2269  CASE DEFAULT
2270  CALL flagerror("Invalid j value",err,error,*999)
2271  END SELECT
2272  CASE DEFAULT
2273  CALL flagerror("Invalid i value",err,error,*999)
2274  END SELECT
2275  ELSE
2276  CALL flagerror("Invalid number of coordinates",err,error,*999)
2277  ENDIF
2278  CASE(coordinate_oblate_spheroidal_type)
2279  CALL flagerror("Not implemented",err,error,*999)
2280  CASE DEFAULT
2281  CALL flagerror("Invalid coordinate type",err,error,*999)
2282  END SELECT
2283 
2284  exits("D2ZX_DP")
2285  RETURN
2286 999 errorsexits("D2ZX_DP",err,error)
2287  RETURN
2288  END FUNCTION d2zx_dp
2289 
2290  !
2291  !================================================================================================================================
2292  !
2293 
2294  !#### Generic-Function: DZX
2295  !### Description:
2296  !### Calculates DZ(:)/DX(J) at X, where Z(:) are rectangalar
2297  !### Cartesian and X(J) are curvilinear coordinates defined by
2298  !### COORDINATE_SYSTEM.
2299  !### Child-functions: DZX_DP,DZX_SP
2300 
2301  !
2302  !================================================================================================================================
2303  !
2304 
2305  FUNCTION dzx_dp(COORDINATE_SYSTEM,I,X,ERR,ERROR)
2306 
2307  !#### Function: DZX_DP
2308  !### Type: REAL(DP)(SIZE(X,1))
2309  !### Description:
2310  !### Calculates DZ(:)/DX(I) at X, where Z(:) are rectangalar
2311  !### Cartesian and X(I) are curvilinear coordinates defined by
2312  !### COORDINATE_SYSTEM.
2313  !### Parent-function: DZX
2314 
2315  !Argument variables
2316  TYPE(coordinate_system_type), INTENT(IN) :: COORDINATE_SYSTEM
2317  INTEGER(INTG), INTENT(IN) :: I
2318  REAL(DP), INTENT(IN) :: X(:)
2319  INTEGER(INTG), INTENT(OUT) :: ERR
2320  TYPE(varying_string), INTENT(OUT) :: ERROR
2321  !Function variable
2322  REAL(DP) :: DZX_DP(size(x,1))
2323  !Local variables
2324  REAL(DP) :: FOCUS
2325 
2326  enters("DZX_DP",err,error,*999)
2327 
2328  dzx_dp=0.0_dp
2329 
2330  IF(SIZE(x,1)<coordinate_system%NUMBER_OF_DIMENSIONS) &
2331  & CALL flagerror("Size of X is less than the number of dimensions", &
2332  & err,error,*999)
2333 
2334  SELECT CASE(coordinate_system%TYPE)
2335  CASE(coordinate_rectangular_cartesian_type)
2336  IF(i>0.AND.i<=coordinate_system%NUMBER_OF_DIMENSIONS) THEN
2337  dzx_dp(i)=1.0_dp
2338  ELSE
2339  CALL flagerror("Invalid i value",err,error,*999)
2340  ENDIF
2341  CASE(coordinate_cylindrical_polar_type)
2342  SELECT CASE(coordinate_system%NUMBER_OF_DIMENSIONS)
2343  CASE(2)
2344  SELECT CASE(i)
2345  CASE(1)
2346  dzx_dp(1)=cos(x(2))
2347  dzx_dp(2)=sin(x(2))
2348  CASE(2)
2349  dzx_dp(1)=-x(1)*sin(x(2))
2350  dzx_dp(2)=x(1)*cos(x(2))
2351  CASE DEFAULT
2352  CALL flagerror("Invalid i value",err,error,*999)
2353  END SELECT
2354  CASE(3)
2355  SELECT CASE(i)
2356  CASE(1)
2357  dzx_dp(1)=cos(x(2))
2358  dzx_dp(2)=sin(x(2))
2359  dzx_dp(3)=0.0_dp
2360  CASE(2)
2361  dzx_dp(1)=-x(1)*sin(x(2))
2362  dzx_dp(2)=x(1)*cos(x(2))
2363  dzx_dp(3)=0.0_dp
2364  CASE(3)
2365  dzx_dp(1)=0.0_dp
2366  dzx_dp(2)=0.0_dp
2367  dzx_dp(3)=1.0_dp
2368  CASE DEFAULT
2369  CALL flagerror("Invalid i value",err,error,*999)
2370  END SELECT
2371  CASE DEFAULT
2372  CALL flagerror("Invalid number of coordinates",err,error,*999)
2373  END SELECT
2374  CASE(coordinate_spherical_polar_type)
2375  IF(coordinate_system%NUMBER_OF_DIMENSIONS==3) THEN
2376  SELECT CASE(i)
2377  CASE(1)
2378  dzx_dp(1)=cos(x(2))*cos(x(3))
2379  dzx_dp(2)=sin(x(2))*cos(x(3))
2380  dzx_dp(3)=sin(x(3))
2381  CASE(2)
2382  dzx_dp(1)=-x(1)*sin(x(2))*cos(x(3))
2383  dzx_dp(2)=x(1)*cos(x(2))*cos(x(3))
2384  dzx_dp(3)=0.0_dp
2385  CASE(3)
2386  dzx_dp(1)=-x(1)*cos(x(2))*sin(x(3))
2387  dzx_dp(2)=-x(1)*sin(x(2))*sin(x(3))
2388  dzx_dp(3)=x(1)*cos(x(3))
2389  CASE DEFAULT
2390  CALL flagerror("Invalid i value",err,error,*999)
2391  END SELECT
2392  ELSE
2393  CALL flagerror("Invalid number of coordinates",err,error,*999)
2394  ENDIF
2395  CASE(coordinate_prolate_spheroidal_type)
2396  IF(coordinate_system%NUMBER_OF_DIMENSIONS==3) THEN
2397  focus=coordinate_system%FOCUS
2398  SELECT CASE(i)
2399  CASE(1)
2400  dzx_dp(1)=focus*sinh(x(1))*cos(x(2))
2401  dzx_dp(2)=focus*cosh(x(1))*sin(x(2))*cos(x(3))
2402  dzx_dp(3)=focus*cosh(x(1))*sin(x(2))*sin(x(3))
2403  CASE(2)
2404  dzx_dp(1)=-focus*cosh(x(1))*sin(x(2))
2405  dzx_dp(2)=focus*sinh(x(1))*cos(x(2))*cos(x(3))
2406  dzx_dp(3)=focus*sinh(x(1))*cos(x(2))*sin(x(3))
2407  CASE(3)
2408  dzx_dp(1)=0.0_dp
2409  dzx_dp(2)=-focus*sinh(x(1))*sin(x(2))*sin(x(3))
2410  dzx_dp(3)=focus*sinh(x(1))*sin(x(2))*cos(x(3))
2411  CASE DEFAULT
2412  CALL flagerror("Invalid i value",err,error,*999)
2413  END SELECT
2414  ELSE
2415  CALL flagerror("Invalid number of coordinates",err,error,*999)
2416  ENDIF
2417  CASE(coordinate_oblate_spheroidal_type)
2418  IF(coordinate_system%NUMBER_OF_DIMENSIONS==3) THEN
2419  focus=coordinate_system%FOCUS
2420  SELECT CASE(i)
2421  CASE(1)
2422  dzx_dp(1)=focus*sinh(x(1))*cos(x(2))*cos(x(3))
2423  dzx_dp(2)=focus*cosh(x(1))*sin(x(2))
2424  dzx_dp(3)=focus*sinh(x(1))*cos(x(2))*sin(x(3))
2425  CASE(2)
2426  dzx_dp(1)=-focus*cosh(x(1))*sin(x(2))*cos(x(3))
2427  dzx_dp(2)=focus*sinh(x(1))*cos(x(2))
2428  dzx_dp(3)=-focus*cosh(x(1))*sin(x(2))*sin(x(3))
2429  CASE(3)
2430  dzx_dp(1)=-focus*cosh(x(1))*cos(x(2))*sin(x(3))
2431  dzx_dp(2)=0.0_dp
2432  dzx_dp(3)=focus*cosh(x(1))*cos(x(2))*cos(x(3))
2433  CASE DEFAULT
2434  CALL flagerror("Invalid i value",err,error,*999)
2435  END SELECT
2436  ELSE
2437  CALL flagerror("Invalid number of coordinates",err,error,*999)
2438  ENDIF
2439  CASE DEFAULT
2440  CALL flagerror("Invalid coordinate type",err,error,*999)
2441  END SELECT
2442 
2443  exits("DZX_DP")
2444  RETURN
2445 999 errorsexits("DZX_DP",err,error)
2446  RETURN
2447  END FUNCTION dzx_dp
2448 
2449  !
2450  !================================================================================================================================
2451  !
2452 
2453  !!TODO:: CHANGE THIS TO A FUNCTION
2454 
2455  !#### Generic-subroutine: COORDINATE_DERIVATIVE_CONVERT_TO_RC
2456  !### Description:
2457  !### COORDINATE_DERIVATIVE_CONVERT_TO_RC performs a coordinate transformation from a
2458  !### coordinate system identified by COORDINATE at the point X
2459  !### with coordinates/derivatives X(nj,nu) to the point with
2460  !### coordinates/derivatives Z(nj) in rectangular cartesian coordinates
2461  !### Child-subroutines: COORDINATE_DERIVATIVE_CONVERT_TO_RC_DP,COORDINATE_DERIVATIVE_CONVERT_TO_RC_SP
2462 
2463  !
2464  !================================================================================================================================
2465  !
2466 
2467  SUBROUTINE coordinate_derivative_convert_to_rc_dp(COORDINATE_SYSTEM,PART_DERIV_TYPE,X,Z,&
2468  & err,error,*)
2469 
2470  !#### Subroutine: COORDINATE_DERIVATIVE_CONVERT_TO_RC_DP
2471  !### Description:
2472  !### COORDINATE_DERIVATIVE_CONVERT_TO_RC_DP performs a coordinate transformation from a
2473  !### coordinate system identified by COORDINATE_SYSTEM at the point X
2474  !### with coordinates/derivatives X(nj,nu) to the point with
2475  !### coordinates/derivatives Z(nj) in rectangular cartesian coordinates
2476  !### for double precision coordinates.
2477  !### Parent-function: COORDINATE_DERIVATIVE_CONVERT_TO_RC
2478 
2479  !Argument variables
2480  TYPE(coordinate_system_type), INTENT(IN) :: COORDINATE_SYSTEM
2481  INTEGER(INTG), INTENT(IN) :: PART_DERIV_TYPE
2482  REAL(DP), INTENT(IN) :: X(:,:)
2483  REAL(DP), INTENT(OUT) :: Z(:)
2484  INTEGER(INTG), INTENT(OUT) :: ERR
2485  TYPE(varying_string), INTENT(OUT) :: ERROR
2486  !Local variables
2487  REAL(DP) :: FOCUS
2488 
2489  enters("COORDINATE_DERIVATIVE_CONVERT_TO_RC_DP",err,error,*999)
2490 
2491 !!TODO: change all second index X(:,?) numbers to their apropriate constant
2492 !!as defined in constants e.g. X(1,2) == X(1,PART_DERIV_S1)
2493 
2494  IF(SIZE(x,1)<coordinate_system%NUMBER_OF_DIMENSIONS) &
2495  & CALL flagerror("Size of X is less than the number of dimensions", &
2496  & err,error,*999)
2497 
2498  IF(SIZE(x,1)==SIZE(z,1)) THEN
2499  SELECT CASE(coordinate_system%TYPE)
2500  CASE(coordinate_rectangular_cartesian_type)
2501  IF(SIZE(x,2)>=part_deriv_type) THEN
2502  z=x(:,part_deriv_type)
2503  ELSE
2504  CALL flagerror("Invalid derivative type",err,error,*999)
2505  ENDIF
2506  CASE(coordinate_cylindrical_polar_type)
2507  SELECT CASE(part_deriv_type)
2508  CASE(no_part_deriv)
2509  z=coordinate_convert_to_rc(coordinate_system,x(:,1),err,error)
2510  IF(err/=0) GOTO 999
2511  CASE(part_deriv_s1)
2512  IF(SIZE(x,2)>=2) THEN
2513  SELECT CASE(coordinate_system%NUMBER_OF_DIMENSIONS)
2514  CASE(2)
2515  z(1)=cos(x(2,1))*x(1,2)-x(1,1)*sin(x(2,1))*x(2,2) !d(x)/d(s1)
2516  z(2)=sin(x(2,1))*x(1,2)+x(1,1)*cos(x(2,1))*x(2,2) !d(y)/d(s1)
2517  CASE(3)
2518  z(1)=cos(x(2,1))*x(1,2)-x(1,1)*sin(x(2,1))*x(2,2) !d(x)/d(s1)
2519  z(2)=sin(x(2,1))*x(1,2)+x(1,1)*cos(x(2,1))*x(2,2) !d(y)/d(s1)
2520  z(3)=x(3,2) !d(z)/d(s1)
2521  CASE DEFAULT
2522  CALL flagerror("Invalid number of coordinates",err,error,*999)
2523  END SELECT
2524  ELSE
2525  CALL flagerror("Not enough X derivatives supplied",err,error,*999)
2526  ENDIF
2527  CASE(part_deriv_s1_s1)
2528  CALL flagerror("Not implemented",err,error,*999)
2529  CASE(part_deriv_s2)
2530  IF(SIZE(x,2)>=4) THEN
2531  SELECT CASE(coordinate_system%NUMBER_OF_DIMENSIONS)
2532  CASE(2)
2533  z(1)=cos(x(2,1))*x(1,4)-x(1,1)*sin(x(2,1))*x(2,4) !d(x)/d(s2)
2534  z(2)=sin(x(2,1))*x(1,4)+x(1,1)*cos(x(2,1))*x(2,4) !d(y)/d(s2)
2535  CASE(3)
2536  z(1)=cos(x(2,1))*x(1,4)-x(1,1)*sin(x(2,1))*x(2,4) !d(x)/d(s2)
2537  z(2)=sin(x(2,1))*x(1,4)+x(1,1)*cos(x(2,1))*x(2,4) !d(y)/d(s2)
2538  z(3)=x(3,4) !d(z)/d(s2)
2539  CASE DEFAULT
2540  CALL flagerror("Invalid number of coordinates",err,error,*999)
2541  END SELECT
2542  ELSE
2543  CALL flagerror("Not enough X derivatives supplied",err,error,*999)
2544  ENDIF
2545  CASE(part_deriv_s2_s2)
2546  CALL flagerror("Not implemented",err,error,*999)
2547  CASE(part_deriv_s1_s2)
2548  IF(SIZE(x,2)>=6) THEN
2549  SELECT CASE(SIZE(x,1))
2550  CASE(2)
2551  z(1)=cos(x(2,1))*x(1,6)-x(1,2)*sin(x(2,1))*x(2,4)-&
2552  & (x(1,4)*sin(x(2,1))*x(2,2)+x(1,1)*cos(x(2,1))*&
2553  & x(2,2)*x(2,4)+x(1,1)*sin(x(2,1))*x(2,6)) !d2(x)/d(s1)d(s2)
2554  z(2)=sin(x(2,1))*x(1,6)+x(1,2)*cos(x(2,1))*x(2,4)+&
2555  & (x(1,4)*cos(x(2,1))*x(2,2)-x(1,1)*sin(x(2,1))*&
2556  & x(2,2)*x(2,4)+x(1,1)*cos(x(2,1))*x(2,6)) !d2(y)/d(s1)d(s2)
2557  CASE(3)
2558  z(1)=cos(x(2,1))*x(1,6)-x(1,2)*sin(x(2,1))*x(2,4)-&
2559  & (x(1,4)*sin(x(2,1))*x(2,2)+x(1,1)*cos(x(2,1))*&
2560  & x(2,2)*x(2,4)+x(1,1)*sin(x(2,1))*x(2,6)) !d2(x)/d(s1)d(s2)
2561  z(2)=sin(x(2,1))*x(1,6)+x(1,2)*cos(x(2,1))*x(2,4)+&
2562  & (x(1,4)*cos(x(2,1))*x(2,2)-x(1,1)*sin(x(2,1))*&
2563  & x(2,2)*x(2,4)+x(1,1)*cos(x(2,1))*x(2,6)) !d2(y)/d(s1)d(s2)
2564  z(3)=x(3,6) !d2(z)/d(s1)d(s2)
2565  CASE DEFAULT
2566  CALL flagerror("Invalid number of coordinates",err,error,*999)
2567  END SELECT
2568  ELSE
2569  CALL flagerror("Not enough X derivatives supplied",err,error,*999)
2570  ENDIF
2571  CASE(part_deriv_s3)
2572  IF(SIZE(x,2)>=7) THEN
2573  SELECT CASE(coordinate_system%NUMBER_OF_DIMENSIONS)
2574  CASE(2)
2575  z(1)=0.0_dp
2576  z(2)=0.0_dp
2577  CASE(3)
2578  z(1)=cos(x(2,1))*x(1,7)-x(1,1)*sin(x(2,1))*x(2,7) !d(x)/d(s3)
2579  z(2)=sin(x(2,1))*x(1,7)+x(1,1)*cos(x(2,1))*x(2,7) !d(y)/d(s3)
2580  z(3)=x(3,7) !d(z)/d(s3)
2581  CASE DEFAULT
2582  CALL flagerror("Invalid number of coordinates",err,error,*999)
2583  END SELECT
2584  ELSE
2585  CALL flagerror("Not enough X derivatives supplied",err,error,*999)
2586  ENDIF
2587  CASE(part_deriv_s3_s3)
2588  CALL flagerror("Not implemented",err,error,*999)
2589  CASE(part_deriv_s1_s3)
2590  IF(SIZE(x,2)>=9) THEN
2591  SELECT CASE(coordinate_system%NUMBER_OF_DIMENSIONS)
2592  CASE(2)
2593  z(1)=0.0_dp
2594  z(2)=0.0_dp
2595  CASE(3)
2596  z(1)=cos(x(2,1))*x(1,9)-x(1,2)*sin(x(2,1))*x(2,7)-&
2597  & (x(1,7)*sin(x(2,1))*x(2,2)+x(1,1)*cos(x(2,1))*&
2598  & x(2,2)*x(2,7)+x(1,1)*sin(x(2,1))*x(2,9)) !d2(x)/d(s1)d(s3)
2599  z(2)=sin(x(2,1))*x(1,9)+x(1,2)*cos(x(2,1))*x(2,7)+&
2600  & (x(1,7)*cos(x(2,1))*x(2,2)-x(1,1)*sin(x(2,1))*&
2601  & x(2,2)*x(2,7)+x(1,1)*cos(x(2,1))*x(2,9)) !d2(y)/d(s1)d(s3)
2602  z(3)=x(3,9) !d2(z)/d(s1)d(s3)
2603  CASE DEFAULT
2604  CALL flagerror("Invalid number of coordinates",err,error,*999)
2605  END SELECT
2606  ELSE
2607  CALL flagerror("Not enough X derivatives supplied",err,error,*999)
2608  ENDIF
2609  CASE(part_deriv_s2_s3)
2610  IF(SIZE(x,2)>=10) THEN
2611  SELECT CASE(coordinate_system%NUMBER_OF_DIMENSIONS)
2612  CASE(2)
2613  z(1)=0.0_dp
2614  z(2)=0.0_dp
2615  CASE(3)
2616  z(1)=cos(x(2,1))*x(1,10)-x(1,4)*sin(x(2,1))*x(2,7)-&
2617  & (x(1,7)*sin(x(2,1))*x(2,4)+x(1,1)*cos(x(2,1))*&
2618  & x(2,4)*x(2,7)+x(1,1)*sin(x(2,1))*x(2,10)) !d2(x)/d(s2)d(s3)
2619  z(2)=sin(x(2,1))*x(1,10)+x(1,4)*cos(x(2,1))*x(2,7)+&
2620  & (x(1,7)*cos(x(2,1))*x(2,4)-x(1,1)*sin(x(2,1))*&
2621  & x(2,4)*x(2,7)+x(1,1)*cos(x(2,1))*x(2,10)) !d2(y)/d(s2)d(s3)
2622  z(3)=x(3,10) !d2(z)/d(s2)d(s3)
2623  CASE DEFAULT
2624  CALL flagerror("Invalid number of coordinates",err,error,*999)
2625  END SELECT
2626  ELSE
2627  CALL flagerror("Not enough X derivatives supplied",err,error,*999)
2628  ENDIF
2629  CASE(part_deriv_s1_s2_s3)
2630  IF(SIZE(x,2)>=11) THEN
2631  SELECT CASE(coordinate_system%NUMBER_OF_DIMENSIONS)
2632  CASE(2)
2633  z(1)=0.0_dp
2634  z(2)=0.0_dp
2635  CASE(3)
2636  z(1)=-cos(x(2,1))*x(2,2)*x(2,4)*x(1,7)-&
2637  & sin(x(2,1))*(x(2,6)*x(1,7)+x(2,4)*x(1,9))-&
2638  & sin(x(2,1))*x(2,2)*x(1,10)+cos(x(2,1))*x(1,11)-&
2639  & cos(x(2,1))*x(2,2)*x(1,4)*x(2,7)-&
2640  & sin(x(2,1))*(x(1,6)*x(2,7)+x(1,4)*x(2,9))+&
2641  & (-cos(x(2,1))*x(1,2)+sin(x(2,1))*x(1,1)*x(2,2))*&
2642  & x(2,4)*x(2,7)-x(1,1)*cos(x(2,1))*(x(2,6)*x(2,7)+&
2643  & x(2,4)*x(2,9))+(-sin(x(2,1))*x(1,2)-x(1,1)*&
2644  & cos(x(2,1))*x(2,2))*x(2,10)-x(1,1)*&
2645  & sin(x(2,1))*x(2,11) !d3(x)/d(s1)d(s2)d(s3)
2646  z(2)=-sin(x(2,1))*x(2,2)*x(2,4)*x(1,7)+&
2647  & cos(x(2,1))*(x(2,6)*x(1,7)+x(2,4)*x(1,9))+&
2648  & cos(x(2,1))*x(2,2)*x(1,10)+sin(x(2,1))*x(1,11)-&
2649  & sin(x(2,1))*x(2,2)*x(1,4)*x(2,7)+&
2650  & cos(x(2,1))*(x(1,6)*x(2,7)+x(1,4)*x(2,9))+&
2651  & (-sin(x(2,1))*x(1,2)-x(1,1)*cos(x(2,1))*x(2,2))*&
2652  & x(2,4)*x(2,7)-x(1,1)*sin(x(2,1))*(x(2,6)*x(2,7)+&
2653  & x(2,4)*x(2,9))+(cos(x(2,1))*x(1,2)-x(1,1)*&
2654  & sin(x(2,1))*x(2,2))*x(2,10)+x(1,1)*&
2655  & cos(x(2,1))*x(2,11) !d3(y)/d(s1)d(s2)d(s3)
2656  z(3)=x(3,11) !d3(z)/d(s1)d(s2)d(s3)
2657  CASE DEFAULT
2658  CALL flagerror("Invalid number of coordinates",err,error,*999)
2659  END SELECT
2660  ELSE
2661  CALL flagerror("Not enough X derivatives supplied",err,error,*999)
2662  ENDIF
2663  CASE DEFAULT
2664  CALL flagerror("Invalid partial derivative type",err,error,*999)
2665  END SELECT
2666  CASE(coordinate_spherical_polar_type)
2667  IF(coordinate_system%NUMBER_OF_DIMENSIONS==3) THEN
2668  SELECT CASE(part_deriv_type)
2669  CASE(no_part_deriv)
2670  z=coordinate_convert_to_rc(coordinate_system,x(:,1),err,error)
2671  IF(err/=0) GOTO 999
2672  CASE(part_deriv_s1,part_deriv_s1_s1,part_deriv_s2,part_deriv_s2_s2,&
2673  & part_deriv_s1_s2,part_deriv_s3,part_deriv_s3_s3,&
2674  & part_deriv_s1_s3,part_deriv_s2_s3,part_deriv_s1_s2_s3)
2675  CALL flagerror("Not implemented",err,error,*999)
2676  CASE DEFAULT
2677  CALL flagerror("Invalid partial derivative type",err,error,*999)
2678  END SELECT
2679  ELSE
2680  CALL flagerror("Invalid number of coordinates",err,error,*999)
2681  ENDIF
2682  CASE(coordinate_prolate_spheroidal_type)
2683  IF(coordinate_system%NUMBER_OF_DIMENSIONS==3) THEN
2684  focus=coordinate_system%FOCUS
2685  SELECT CASE(part_deriv_type)
2686  CASE(no_part_deriv)
2687  z=coordinate_convert_to_rc(coordinate_system,x(:,1),err,error)
2688  IF(err/=0) GOTO 999
2689  CASE(part_deriv_s1)
2690  IF(SIZE(x,2)>=2) THEN
2691  z(1)=focus*sinh(x(1,1))*cos(x(2,1))*x(1,2)-&
2692  & focus*cosh(x(1,1))*sin(x(2,1))*x(2,2) !d(x)/d(s1)
2693  z(2)=focus*cosh(x(1,1))*sin(x(2,1))*cos(x(3,1))*x(1,2)+&
2694  & focus*sinh(x(1,1))*cos(x(2,1))*cos(x(3,1))*x(2,2)-&
2695  & focus*sinh(x(1,1))*sin(x(2,1))*sin(x(3,1))*&
2696  & x(3,2) !d(y)/d(s1)
2697  z(3)=focus*cosh(x(1,1))*sin(x(2,1))*sin(x(3,1))*x(1,2)+&
2698  & focus*sinh(x(1,1))*cos(x(2,1))*sin(x(3,1))*x(2,2)+&
2699  & focus*sinh(x(1,1))*sin(x(2,1))*cos(x(3,1))*&
2700  & x(3,2) !d(z)/d(s1)
2701  ELSE
2702  CALL flagerror("Invalid derivative type",err,error,*999)
2703  ENDIF
2704  CASE(part_deriv_s1_s1)
2705  CALL flagerror("Not implemented",err,error,*999)
2706  CASE(part_deriv_s2)
2707  IF(SIZE(x,2)>=4) THEN
2708  z(1)=focus*sinh(x(1,1))*cos(x(2,1))*x(1,4)-&
2709  & focus*cosh(x(1,1))*sin(x(2,1))*x(2,4) !d(x)/d(s2)
2710  z(2)=focus*cosh(x(1,1))*sin(x(2,1))*cos(x(3,1))*x(1,4)+&
2711  & focus*sinh(x(1,1))*cos(x(2,1))*cos(x(3,1))*x(2,4)-&
2712  & focus*sinh(x(1,1))*sin(x(2,1))*sin(x(3,1))*&
2713  & x(3,4) !d(y)/d(s2)
2714  z(3)=focus*cosh(x(1,1))*sin(x(2,1))*sin(x(3,1))*x(1,4)+&
2715  & focus*sinh(x(1,1))*cos(x(2,1))*sin(x(3,1))*x(2,4)+&
2716  & focus*sinh(x(1,1))*sin(x(2,1))*cos(x(3,1))*&
2717  & x(3,4) !d(z)/d(s2)
2718  ELSE
2719  CALL flagerror("Not enough X derivatives supplied",&
2720  & err,error,*999)
2721  ENDIF
2722  CASE(part_deriv_s2_s2)
2723  CALL flagerror("Not implemented",err,error,*999)
2724  CASE(part_deriv_s1_s2)
2725  IF(SIZE(x,2)>=6) THEN
2726  z(1)=focus*(sinh(x(1,1))*cos(x(2,1))*x(1,6)+&
2727  & x(1,2)*(cosh(x(1,1))*cos(x(2,1))*x(1,4)-&
2728  & sinh(x(1,1))*sin(x(2,1))*x(2,4))-&
2729  & cosh(x(1,1))*sin(x(2,1))*x(2,6)-&
2730  & x(2,2)*(sinh(x(1,1))*sin(x(2,1))*x(1,4)+&
2731  & cosh(x(1,1))*cos(x(2,1))*x(2,4))) !d2(x)/d(s1)d(s2)
2732  z(2)=focus*(cosh(x(1,1))*sin(x(2,1))*cos(x(3,1))*x(1,6)+&
2733  & x(1,2)*(sinh(x(1,1))*sin(x(2,1))*cos(x(3,1))*&
2734  & x(1,4)+cosh(x(1,1))*cos(x(2,1))*cos(x(3,1))*&
2735  & x(2,4)-cosh(x(1,1))*sin(x(2,1))*sin(x(3,1))*&
2736  & x(3,4))+sinh(x(1,1))*cos(x(2,1))*cos(x(3,1))*&
2737  & x(2,6)+x(2,2)*(cosh(x(1,1))*cos(x(2,1))*cos(x(3,1))*&
2738  & x(1,4)-sinh(x(1,1))*sin(x(2,1))*cos(x(3,1))*x(2,4)-&
2739  & sinh(x(1,1))*cos(x(2,1))*sin(x(3,1))*x(3,4))-&
2740  & sinh(x(1,1))*sin(x(2,1))*sin(x(3,1))*x(3,6)-&
2741  & x(3,2)*(cosh(x(1,1))*sin(x(2,1))*sin(x(3,1))*&
2742  & x(1,4)+sinh(x(1,1))*cos(x(2,1))*sin(x(3,1))*&
2743  & x(2,4)+sinh(x(1,1))*sin(x(2,1))*cos(x(3,1))*&
2744  & x(3,4))) !d2(y)/d(s1)d(s2)
2745  z(3)=focus*(cosh(x(1,1))*sin(x(2,1))*sin(x(3,1))*x(1,6)+&
2746  & x(1,2)*(sinh(x(1,1))*sin(x(2,1))*sin(x(3,1))*&
2747  & x(1,4)+cosh(x(1,1))*cos(x(2,1))*sin(x(3,1))*&
2748  & x(2,4)+cosh(x(1,1))*sin(x(2,1))*cos(x(3,1))*&
2749  & x(3,4))+sinh(x(1,1))*cos(x(2,1))*sin(x(3,1))*x(2,6)+&
2750  & x(2,2)*(cosh(x(1,1))*cos(x(2,1))*sin(x(3,1))*&
2751  & x(1,4)-sinh(x(1,1))*sin(x(2,1))*sin(x(3,1))*&
2752  & x(2,4)+sinh(x(1,1))*cos(x(2,1))*cos(x(3,1))*&
2753  & x(3,4))+sinh(x(1,1))*sin(x(2,1))*cos(x(3,1))*x(3,6)+&
2754  & x(3,2)*(cosh(x(1,1))*sin(x(2,1))*cos(x(3,1))*&
2755  & x(1,4)+sinh(x(1,1))*cos(x(2,1))*cos(x(3,1))*&
2756  & x(2,4)-sinh(x(1,1))*sin(x(2,1))*sin(x(3,1))*&
2757  & x(3,4))) !d2(z)/d(s1)d(s2)
2758  ELSE
2759  CALL flagerror("Not enough X derivatives supplied",&
2760  & err,error,*999)
2761  ENDIF
2762  CASE(part_deriv_s3)
2763  IF(SIZE(x,2)>=7) THEN
2764  z(1)=focus*sinh(x(1,1))*cos(x(2,1))*x(1,7)-&
2765  & focus*cosh(x(1,1))*sin(x(2,1))*x(2,7) !d(x)/d(s3)
2766  z(2)=focus*cosh(x(1,1))*sin(x(2,1))*cos(x(3,1))*x(1,7)+&
2767  & focus*sinh(x(1,1))*cos(x(2,1))*cos(x(3,1))*x(2,7)-&
2768  & focus*sinh(x(1,1))*sin(x(2,1))*sin(x(3,1))*&
2769  & x(3,7) !d(y)/d(s3)
2770  z(3)=focus*cosh(x(1,1))*sin(x(2,1))*sin(x(3,1))*x(1,7)+&
2771  & focus*sinh(x(1,1))*cos(x(2,1))*sin(x(3,1))*x(2,7)+&
2772  & focus*sinh(x(1,1))*sin(x(2,1))*cos(x(3,1))*&
2773  & x(3,7) !d(z)/d(s3)
2774  ELSE
2775  CALL flagerror("Not enough X derivatives supplied",&
2776  & err,error,*999)
2777  ENDIF
2778  CASE(part_deriv_s3_s3)
2779  CALL flagerror("Not implemented",err,error,*999)
2780  CASE(part_deriv_s1_s3)
2781  IF(SIZE(x,2)>=9) THEN
2782  z(1)=focus*(sinh(x(1,1))*cos(x(2,1))*x(1,9)+&
2783  & x(1,2)*(cosh(x(1,1))*cos(x(2,1))*x(1,7)-&
2784  & sinh(x(1,1))*sin(x(2,1))*x(2,7))-&
2785  & cosh(x(1,1))*sin(x(2,1))*x(2,9)-&
2786  & x(2,2)*(sinh(x(1,1))*sin(x(2,1))*x(1,7)+&
2787  & cosh(x(1,1))*cos(x(2,1))*x(2,7))) !d2(x)/d(s1)d(s3)
2788  z(2)=focus*(cosh(x(1,1))*sin(x(2,1))*cos(x(3,1))*x(1,9)+&
2789  & x(1,2)*(sinh(x(1,1))*sin(x(2,1))*cos(x(3,1))*x(1,7)+&
2790  & cosh(x(1,1))*cos(x(2,1))*cos(x(3,1))*x(2,7)-&
2791  & cosh(x(1,1))*sin(x(2,1))*sin(x(3,1))*x(3,7))+&
2792  & sinh(x(1,1))*cos(x(2,1))*cos(x(3,1))*x(2,9)+&
2793  & x(2,2)*(cosh(x(1,1))*cos(x(2,1))*cos(x(3,1))*x(1,7)-&
2794  & sinh(x(1,1))*sin(x(2,1))*cos(x(3,1))*x(2,7)-&
2795  & sinh(x(1,1))*cos(x(2,1))*sin(x(3,1))*x(3,7))-&
2796  & sinh(x(1,1))*sin(x(2,1))*sin(x(3,1))*x(3,9)-&
2797  & x(3,2)*(cosh(x(1,1))*sin(x(2,1))*sin(x(3,1))*x(1,7)+&
2798  & sinh(x(1,1))*cos(x(2,1))*sin(x(3,1))* x(2,7)+&
2799  & sinh(x(1,1))*sin(x(2,1))*cos(x(3,1))*&
2800  & x(3,7))) !d2(y)/d(s1)d(s3)
2801  z(3)=focus*(cosh(x(1,1))*sin(x(2,1))*sin(x(3,1))*x(1,9)+&
2802  & x(1,2)*(sinh(x(1,1))*sin(x(2,1))*sin(x(3,1))*x(1,7)+&
2803  & cosh(x(1,1))*cos(x(2,1))*sin(x(3,1))*x(2,7)+&
2804  & cosh(x(1,1))*sin(x(2,1))*cos(x(3,1))*x(3,7))+&
2805  & sinh(x(1,1))*cos(x(2,1))*sin(x(3,1))*x(2,9)+&
2806  & x(2,2)*(cosh(x(1,1))*cos(x(2,1))*sin(x(3,1))*x(1,7)-&
2807  & sinh(x(1,1))*sin(x(2,1))*sin(x(3,1))*x(2,7)+&
2808  & sinh(x(1,1))*cos(x(2,1))*cos(x(3,1))*x(3,7))+&
2809  & sinh(x(1,1))*sin(x(2,1))*cos(x(3,1))*x(3,9)+&
2810  & x(3,2)*(cosh(x(1,1))*sin(x(2,1))*cos(x(3,1))*x(1,7)+&
2811  & sinh(x(1,1))*cos(x(2,1))*cos(x(3,1))*x(2,7)-&
2812  & sinh(x(1,1))*sin(x(2,1))*sin(x(3,1))*&
2813  & x(3,7))) !d2(z)/d(s1)d(s3)
2814  ELSE
2815  CALL flagerror("Not enough X derivatives supplied",&
2816  & err,error,*999)
2817  ENDIF
2818  CASE(part_deriv_s2_s3)
2819  IF(SIZE(x,2)>=10) THEN
2820  z(1)=focus*(sinh(x(1,1))*cos(x(2,1))*x(1,10)+&
2821  & x(1,4)*(cosh(x(1,1))*cos(x(2,1))*x(1,7)-&
2822  & sinh(x(1,1))*sin(x(2,1))*x(2,7))-&
2823  & cosh(x(1,1))*sin(x(2,1))*x(2,10)-&
2824  & x(2,4)*(sinh(x(1,1))*sin(x(2,1))*x(1,7)+&
2825  & cosh(x(1,1))*cos(x(2,1))*x(2,7))) !d2(x)/d(s2)d(s3)
2826  z(2)=focus*(cosh(x(1,1))*sin(x(2,1))*cos(x(3,1))*x(1,10)+&
2827  & x(1,4)*(sinh(x(1,1))*sin(x(2,1))*cos(x(3,1))*x(1,7)+&
2828  & cosh(x(1,1))*cos(x(2,1))*cos(x(3,1))*x(2,7)-&
2829  & cosh(x(1,1))*sin(x(2,1))*sin(x(3,1))*x(3,7))+&
2830  & sinh(x(1,1))*cos(x(2,1))*cos(x(3,1))*x(2,10)+&
2831  & x(2,4)*(cosh(x(1,1))*cos(x(2,1))*cos(x(3,1))*x(1,7)-&
2832  & sinh(x(1,1))*sin(x(2,1))*cos(x(3,1))*x(2,7)-&
2833  & sinh(x(1,1))*cos(x(2,1))*sin(x(3,1))*x(3,7))-&
2834  & sinh(x(1,1))*sin(x(2,1))*sin(x(3,1))*x(3,10)-&
2835  & x(3,4)*(cosh(x(1,1))*sin(x(2,1))*sin(x(3,1))*x(1,7)+&
2836  & sinh(x(1,1))*cos(x(2,1))*sin(x(3,1))*x(2,7)+&
2837  & sinh(x(1,1))*sin(x(2,1))*cos(x(3,1))*&
2838  & x(3,7))) !d2(y)/d(s2)d(s3)
2839  z(3)=focus*(cosh(x(1,1))*sin(x(2,1))*sin(x(3,1))*x(1,10)+&
2840  & x(1,4)*(sinh(x(1,1))*sin(x(2,1))*sin(x(3,1))*x(1,7)+&
2841  & cosh(x(1,1))*cos(x(2,1))*sin(x(3,1))*x(2,7)+&
2842  & cosh(x(1,1))*sin(x(2,1))*cos(x(3,1))*x(3,7))+&
2843  & sinh(x(1,1))*cos(x(2,1))*sin(x(3,1))*x(2,10)+&
2844  & x(2,4)*(cosh(x(1,1))*cos(x(2,1))*sin(x(3,1))*x(1,7)-&
2845  & sinh(x(1,1))*sin(x(2,1))*sin(x(3,1))*x(2,7)+&
2846  & sinh(x(1,1))*cos(x(2,1))*cos(x(3,1))*x(3,7))+&
2847  & sinh(x(1,1))*sin(x(2,1))*cos(x(3,1))*x(3,10)+&
2848  & x(3,4)*(cosh(x(1,1))*sin(x(2,1))*cos(x(3,1))*x(1,7)+&
2849  & sinh(x(1,1))*cos(x(2,1))*cos(x(3,1))*x(2,7)-&
2850  & sinh(x(1,1))*sin(x(2,1))*sin(x(3,1))*&
2851  & x(3,7))) !d2(z)/d(s2)d(s3)
2852  ELSE
2853  CALL flagerror("Not enough X derivatives supplied",&
2854  & err,error,*999)
2855  ENDIF
2856  CASE(part_deriv_s1_s2_s3)
2857  IF(SIZE(x,2)>=11) THEN
2858  z(1)=focus*((sinh(x(1,1))*cos(x(2,1))*x(1,2)-&
2859  & cosh(x(1,1))*sin(x(2,1))*x(2,2))*x(1,4)*x(1,7)+&
2860  & cosh(x(1,1))*cos(x(2,1))*(x(1,6)*x(1,7)+x(1,4)*x(1,9))+&
2861  & (-cosh(x(1,1))*sin(x(2,1))*x(1,2)-&
2862  & sinh(x(1,1))*cos(x(2,1))*x(2,2))*x(2,4)*x(1,7)-&
2863  & sinh(x(1,1))*sin(x(2,1))*(x(2,6)*x(1,7)+x(2,4)*x(1,9))+&
2864  & (cosh(x(1,1))*cos(x(2,1))*x(1,2)-&
2865  & sinh(x(1,1))*sin(x(2,1))*x(2,2))*x(1,10)+&
2866  & sinh(x(1,1))*cos(x(2,1))*x(1,11)+&
2867  & (-cosh(x(1,1))*sin(x(2,1))*x(1,2)-&
2868  & sinh(x(1,1))*cos(x(2,1))*x(2,2))*x(1,4)*x(2,7)-&
2869  & sinh(x(1,1))*sin(x(2,1))*(x(1,6)*x(2,7)+x(1,4)*x(2,9))+&
2870  & (-sinh(x(1,1))*cos(x(2,1))*x(1,2)+&
2871  & cosh(x(1,1))*sin(x(2,1))*x(2,2))*x(2,4)*x(2,7)-&
2872  & cosh(x(1,1))*cos(x(2,1))*(x(2,6)*x(2,7)+x(2,4)*x(2,9))+&
2873  & (-sinh(x(1,1))*sin(x(2,1))*x(1,2)-&
2874  & cosh(x(1,1))*cos(x(2,1))*x(2,2))*x(2,10)-&
2875  & cosh(x(1,1))*sin(x(2,1))*x(2,11)) !d3(x)/d(s1)d(s2)d(s3)
2876  z(2)=focus*((cosh(x(1,1))*sin(x(2,1))*cos(x(3,1))*x(1,2)+&
2877  & sinh(x(1,1))*cos(x(2,1))*cos(x(3,1))*x(2,2)-&
2878  & sinh(x(1,1))*sin(x(2,1))*sin(x(3,1))*x(3,2))*&
2879  & x(1,4)*x(1,7)+sinh(x(1,1))*sin(x(2,1))*cos(x(3,1))*&
2880  & (x(1,6)*x(1,7)+x(1,4)*x(1,9))+&
2881  & (sinh(x(1,1))*cos(x(2,1))*cos(x(3,1))*x(1,2)-&
2882  & cosh(x(1,1))*sin(x(2,1))*cos(x(3,1))*x(2,2)-&
2883  & cosh(x(1,1))*cos(x(2,1))*sin(x(3,1))*x(3,2))*&
2884  & x(2,4)*x(1,7)+cosh(x(1,1))*cos(x(2,1))*cos(x(3,1))*&
2885  & (x(2,6)*x(1,7)+x(2,4)*x(1,9))+&
2886  & (-sinh(x(1,1))*sin(x(2,1))*sin(x(3,1))*x(1,2)-&
2887  & cosh(x(1,1))*cos(x(2,1))*sin(x(3,1))*x(2,2)-&
2888  & cosh(x(1,1))*sin(x(2,1))*cos(x(3,1))*x(3,2))*&
2889  & x(3,4)*x(1,7)-cosh(x(1,1))*sin(x(2,1))*sin(x(3,1))*&
2890  & (x(3,6)*x(1,7)+x(3,4)*x(1,9))+&
2891  & (sinh(x(1,1))*sin(x(2,1))*cos(x(3,1))*x(1,2)+&
2892  & cosh(x(1,1))*cos(x(2,1))*cos(x(3,1))*x(2,2)-&
2893  & cosh(x(1,1))*sin(x(2,1))*sin(x(3,1))*x(3,2))*x(1,10)+&
2894  & cosh(x(1,1))*sin(x(2,1))*cos(x(3,1))*x(1,11))+&
2895  & focus*((sinh(x(1,1))*cos(x(2,1))*cos(x(3,1))*x(1,2)-&
2896  & cosh(x(1,1))*sin(x(2,1))*cos(x(3,1))*x(2,2)-&
2897  & cosh(x(1,1))*cos(x(2,1))*sin(x(3,1))*x(3,2))*&
2898  & x(1,4)*x(2,7)+cosh(x(1,1))*cos(x(2,1))*cos(x(3,1))*&
2899  & (x(1,6)*x(2,7)+x(1,4)*x(2,9))+&
2900  & (-cosh(x(1,1))*sin(x(2,1))*cos(x(3,1))*x(1,2)-&
2901  & sinh(x(1,1))*cos(x(2,1))*cos(x(3,1))*x(2,2)+&
2902  & sinh(x(1,1))*sin(x(2,1))*sin(x(3,1))*x(3,2))*&
2903  & x(2,4)*x(2,7)-sinh(x(1,1))*sin(x(2,1))*cos(x(3,1))*&
2904  & (x(2,6)*x(2,7)+x(2,4)*x(2,9))+&
2905  & (-cosh(x(1,1))*cos(x(2,1))*sin(x(3,1))*x(1,2)+&
2906  & sinh(x(1,1))*sin(x(2,1))*sin(x(3,1))*x(2,2)-&
2907  & sinh(x(1,1))*cos(x(2,1))*cos(x(3,1))*x(3,2))*&
2908  & x(3,4)*x(2,7)-sinh(x(1,1))*cos(x(2,1))*sin(x(3,1))*&
2909  & (x(3,6)*x(2,7)+x(3,4)*x(2,9))+&
2910  & (cosh(x(1,1))*cos(x(2,1))*cos(x(3,1))*x(1,2)-&
2911  & sinh(x(1,1))*sin(x(2,1))*cos(x(3,1))*x(2,2)-&
2912  & sinh(x(1,1))*cos(x(2,1))*sin(x(3,1))*x(3,2))*x(2,10)+&
2913  & sinh(x(1,1))*cos(x(2,1))*cos(x(3,1))*x(2,11))+&
2914  & focus*((-sinh(x(1,1))*sin(x(2,1))*sin(x(3,1))*x(1,2)-&
2915  & cosh(x(1,1))*cos(x(2,1))*sin(x(3,1))*x(2,2)-&
2916  & cosh(x(1,1))*sin(x(2,1))*cos(x(3,1))*x(3,2))*&
2917  & x(1,4)*x(3,7)-cosh(x(1,1))*sin(x(2,1))*sin(x(3,1))*&
2918  & (x(1,6)*x(3,7)+x(1,4)*x(3,9))+&
2919  & (-cosh(x(1,1))*cos(x(2,1))*sin(x(3,1))*x(1,2)+&
2920  & sinh(x(1,1))*sin(x(2,1))*sin(x(3,1))*x(2,2)-&
2921  & sinh(x(1,1))*cos(x(2,1))*cos(x(3,1))*x(3,2))*&
2922  & x(2,4)*x(3,7)-sinh(x(1,1))*cos(x(2,1))*sin(x(3,1))*&
2923  & (x(2,6)*x(3,7)+x(2,4)*x(3,9))+&
2924  & (-cosh(x(1,1))*sin(x(2,1))*cos(x(3,1))*x(1,2)-&
2925  & sinh(x(1,1))*cos(x(2,1))*cos(x(3,1))*x(2,2)+&
2926  & sinh(x(1,1))*sin(x(2,1))*sin(x(3,1))*x(3,2))*&
2927  & x(3,4)*x(3,7)-sinh(x(1,1))*sin(x(2,1))*cos(x(3,1))*&
2928  & (x(3,6)*x(3,7)+x(3,4)*x(3,9))+&
2929  & (-cosh(x(1,1))*sin(x(2,1))*sin(x(3,1))*x(1,2)-&
2930  & sinh(x(1,1))*cos(x(2,1))*sin(x(3,1))*x(2,2)-&
2931  & sinh(x(1,1))*sin(x(2,1))*cos(x(3,1))*x(3,2))*x(3,10)-&
2932  & sinh(x(1,1))*sin(x(2,1))*sin(x(3,1))*&
2933  & x(3,11)) !d3(y)/d(s1)d(s2)d(s3)
2934  z(3)=focus*((cosh(x(1,1))*sin(x(2,1))*sin(x(3,1))*x(1,2)+&
2935  & sinh(x(1,1))*cos(x(2,1))*sin(x(3,1))*x(2,2)+&
2936  & sinh(x(1,1))*sin(x(2,1))*cos(x(3,1))*x(3,2))*&
2937  & x(1,4)*x(1,7)+sinh(x(1,1))*sin(x(2,1))*sin(x(3,1))*&
2938  & (x(1,6)*x(1,7)+x(1,4)*x(1,9))+&
2939  & (sinh(x(1,1))*cos(x(2,1))*sin(x(3,1))*x(1,2)-&
2940  & cosh(x(1,1))*sin(x(2,1))*sin(x(3,1))*x(2,2)+&
2941  & cosh(x(1,1))*cos(x(2,1))*cos(x(3,1))*x(3,2))*&
2942  & x(2,4)*x(1,7)+cosh(x(1,1))*cos(x(2,1))*sin(x(3,1))*&
2943  & (x(2,6)*x(1,7)+x(2,4)*x(1,9))+&
2944  & (sinh(x(1,1))*sin(x(2,1))*cos(x(3,1))*x(1,2)+&
2945  & cosh(x(1,1))*cos(x(2,1))*cos(x(3,1))*x(2,2)-&
2946  & cosh(x(1,1))*sin(x(2,1))*sin(x(3,1))*x(3,2))*&
2947  & x(3,4)*x(1,7)+cosh(x(1,1))*sin(x(2,1))*cos(x(3,1))*&
2948  & (x(3,6)*x(1,7)+x(3,4)*x(1,9))+&
2949  & (sinh(x(1,1))*sin(x(2,1))*sin(x(3,1))*x(1,2)+&
2950  & cosh(x(1,1))*cos(x(2,1))*sin(x(3,1))*x(2,2)+&
2951  & cosh(x(1,1))*sin(x(2,1))*cos(x(3,1))*x(3,2))*x(1,10)+&
2952  & cosh(x(1,1))*sin(x(2,1))*sin(x(3,1))*x(1,11))+&
2953  & focus*((sinh(x(1,1))*cos(x(2,1))*sin(x(3,1))*x(1,2)-&
2954  & cosh(x(1,1))*sin(x(2,1))*sin(x(3,1))*x(2,2)+&
2955  & cosh(x(1,1))*cos(x(2,1))*cos(x(3,1))*x(3,2))*&
2956  & x(1,4)*x(2,7)+cosh(x(1,1))*cos(x(2,1))*sin(x(3,1))*&
2957  & (x(1,6)*x(2,7)+x(1,4)*x(2,9))+&
2958  & (-cosh(x(1,1))*sin(x(2,1))*sin(x(3,1))*x(1,2)-&
2959  & sinh(x(1,1))*cos(x(2,1))*sin(x(3,1))*x(2,2)-&
2960  & sinh(x(1,1))*sin(x(2,1))*cos(x(3,1))*x(3,2))*&
2961  & x(2,4)*x(2,7)-sinh(x(1,1))*sin(x(2,1))*sin(x(3,1))*&
2962  & (x(2,6)*x(2,7)+x(2,4)*x(2,9))+&
2963  & (cosh(x(1,1))*cos(x(2,1))*cos(x(3,1))*x(1,2)-&
2964  & sinh(x(1,1))*sin(x(2,1))*cos(x(3,1))*x(2,2)-&
2965  & sinh(x(1,1))*cos(x(2,1))*sin(x(3,1))*x(3,2))*&
2966  & x(3,4)*x(2,7)+sinh(x(1,1))*cos(x(2,1))*cos(x(3,1))*&
2967  & (x(3,6)*x(2,7)+x(3,4)*x(2,9))+&
2968  & (cosh(x(1,1))*cos(x(2,1))*sin(x(3,1))*x(1,2)-&
2969  & sinh(x(1,1))*sin(x(2,1))*sin(x(3,1))*x(2,2)+&
2970  & sinh(x(1,1))*cos(x(2,1))*cos(x(3,1))*x(3,2))*x(2,10)+&
2971  & sinh(x(1,1))*cos(x(2,1))*sin(x(3,1))*x(2,11))+&
2972  & focus*((sinh(x(1,1))*sin(x(2,1))*cos(x(3,1))*x(1,2)+&
2973  & cosh(x(1,1))*cos(x(2,1))*cos(x(3,1))*x(2,2)-&
2974  & cosh(x(1,1))*sin(x(2,1))*sin(x(3,1))*x(3,2))*&
2975  & x(1,4)*x(3,7)+cosh(x(1,1))*sin(x(2,1))*cos(x(3,1))*&
2976  & (x(1,6)*x(3,7)+x(1,4)*x(3,9))+&
2977  & (cosh(x(1,1))*cos(x(2,1))*cos(x(3,1))*x(1,2)-&
2978  & sinh(x(1,1))*sin(x(2,1))*cos(x(3,1))*x(2,2)-&
2979  & sinh(x(1,1))*cos(x(2,1))*sin(x(3,1))*x(3,2))*&
2980  & x(2,4)*x(3,7)+sinh(x(1,1))*cos(x(2,1))*cos(x(3,1))*&
2981  & (x(2,6)*x(3,7)+x(2,4)*x(3,9))+&
2982  & (-cosh(x(1,1))*sin(x(2,1))*sin(x(3,1))*x(1,2)-&
2983  & sinh(x(1,1))*cos(x(2,1))*sin(x(3,1))*x(2,2)-&
2984  & sinh(x(1,1))*sin(x(2,1))*cos(x(3,1))*x(3,2))*&
2985  & x(3,4)*x(3,7)-sinh(x(1,1))*sin(x(2,1))*sin(x(3,1))*&
2986  & (x(3,6)*x(3,7)+x(3,4)*x(3,9))+&
2987  & (cosh(x(1,1))*sin(x(2,1))*cos(x(3,1))*x(1,2)+&
2988  & sinh(x(1,1))*cos(x(2,1))*cos(x(3,1))*x(2,2)-&
2989  & sinh(x(1,1))*sin(x(2,1))*sin(x(3,1))*x(3,2))*x(3,10)+&
2990  & sinh(x(1,1))*sin(x(2,1))*cos(x(3,1))*&
2991  & x(3,11)) !d3(z)/d(s1)d(s2)d(s3)
2992  ELSE
2993  CALL flagerror("Not enough X derivatives supplied",&
2994  & err,error,*999)
2995  ENDIF
2996  CASE DEFAULT
2997  CALL flagerror("Invalid partial derivative type",err,error,*999)
2998  END SELECT
2999  ENDIF
3000  CASE(coordinate_oblate_spheroidal_type)
3001  IF(coordinate_system%NUMBER_OF_DIMENSIONS==3) THEN
3002  focus=coordinate_system%FOCUS
3003  SELECT CASE(part_deriv_type)
3004  CASE(no_part_deriv)
3005  z=coordinate_convert_to_rc(coordinate_system,x(:,1),err,error)
3006  IF(err/=0) GOTO 999
3007  CASE(part_deriv_s1,part_deriv_s1_s1,part_deriv_s2,part_deriv_s2_s2,&
3008  & part_deriv_s1_s2,part_deriv_s3,part_deriv_s3_s3,&
3009  & part_deriv_s1_s3,part_deriv_s2_s3,part_deriv_s1_s2_s3)
3010  CALL flagerror("Not implemented",err,error,*999)
3011  CASE DEFAULT
3012  CALL flagerror("Invalid partial derivative type",err,error,*999)
3013  END SELECT
3014  ELSE
3015  CALL flagerror("Invalid number of coordinates",err,error,*999)
3016  ENDIF
3017  CASE DEFAULT
3018  CALL flagerror("Invalid coordinate type",err,error,*999)
3019  END SELECT
3020  ELSE
3021  CALL flagerror("The sizes of the vectors X and Z do not match",&
3022  & err,error,*999)
3023  ENDIF
3024 
3025  exits("COORDINATE_DERIVATIVE_CONVERT_TO_RC_DP")
3026  RETURN
3027 999 errorsexits("COORDINATE_DERIVATIVE_CONVERT_TO_RC_DP",err,error)
3028  RETURN 1
3029  END SUBROUTINE coordinate_derivative_convert_to_rc_dp
3030 
3031  !
3032  !================================================================================================================================
3033  !
3034 
3035  SUBROUTINE coordinate_derivative_convert_to_rc_sp(COORDINATE_SYSTEM,PART_DERIV_TYPE,X,Z,&
3036  & err,error,*)
3037 
3038  !#### Subroutine: COORDINATE_DERIVATIVE_CONVERT_TO_RC_SP
3039  !### Description:
3040  !### COORDINATE_DERIVATIVE_CONVERT_TO_RC_SP performs a coordinate transformation from a
3041  !### coordinate system identified by COORDINATE_SYSTEM at the point X
3042  !### with coordinates/derivatives X(nj,nu) to the point with
3043  !### coordinates/derivatives Z(nj) in rectangular cartesian coordinates
3044  !### for single precision coordinates.
3045  !### Parent-function: COORDINATE_DERIVATIVE_CONVERT_TO_RC
3046 
3047  !Argument variables
3048  TYPE(coordinate_system_type), INTENT(IN) :: COORDINATE_SYSTEM
3049  INTEGER(INTG), INTENT(IN) :: PART_DERIV_TYPE
3050  REAL(SP), INTENT(IN) :: X(:,:)
3051  REAL(SP), INTENT(OUT) :: Z(:)
3052  INTEGER(INTG), INTENT(OUT) :: ERR
3053  TYPE(varying_string), INTENT(OUT) :: ERROR
3054  !Local variables
3055  REAL(SP) :: FOCUS
3056 
3057  enters("COORDINATE_DERIVATIVE_CONVERT_TO_RC_SP",err,error,*999)
3058 
3059 !!TODO: change all second index X(:,?) numbers to their apropriate constant
3060 !!as defined in constants e.g. X(1,2) == X(1,PART_DERIV_S1)
3061 
3062  IF(SIZE(x,1)<coordinate_system%NUMBER_OF_DIMENSIONS) &
3063  & CALL flagerror("Size of X is less than the number of dimensions", &
3064  & err,error,*999)
3065 
3066  IF(SIZE(x,1)==SIZE(z,1)) THEN
3067  SELECT CASE(coordinate_system%TYPE)
3068  CASE(coordinate_rectangular_cartesian_type)
3069  IF(SIZE(x,2)>=part_deriv_type) THEN
3070  z=x(:,part_deriv_type)
3071  ELSE
3072  CALL flagerror("Invalid partial derivative type",err,error,*999)
3073  ENDIF
3074  CASE(coordinate_cylindrical_polar_type)
3075  SELECT CASE(part_deriv_type)
3076  CASE(no_part_deriv)
3077  z=coordinate_convert_to_rc(coordinate_system,x(:,1),err,error)
3078  IF(err/=0) GOTO 999
3079  CASE(part_deriv_s1)
3080  IF(SIZE(x,2)>=2) THEN
3081  SELECT CASE(coordinate_system%NUMBER_OF_DIMENSIONS)
3082  CASE(2)
3083  z(1)=cos(x(2,1))*x(1,2)-x(1,1)*sin(x(2,1))*x(2,2) !d(x)/d(s1)
3084  z(2)=sin(x(2,1))*x(1,2)+x(1,1)*cos(x(2,1))*x(2,2) !d(y)/d(s1)
3085  CASE(3)
3086  z(1)=cos(x(2,1))*x(1,2)-x(1,1)*sin(x(2,1))*x(2,2) !d(x)/d(s1)
3087  z(2)=sin(x(2,1))*x(1,2)+x(1,1)*cos(x(2,1))*x(2,2) !d(y)/d(s1)
3088  z(3)=x(3,2) !d(z)/d(s1)
3089  CASE DEFAULT
3090  CALL flagerror("Invalid number of coordinates",err,error,*999)
3091  END SELECT
3092  ELSE
3093  CALL flagerror("Not enough X derivatives supplied",err,error,*999)
3094  ENDIF
3095  CASE(part_deriv_s1_s1)
3096  CALL flagerror("Not implemented",err,error,*999)
3097  CASE(part_deriv_s2)
3098  IF(SIZE(x,2)>=4) THEN
3099  SELECT CASE(coordinate_system%NUMBER_OF_DIMENSIONS)
3100  CASE(2)
3101  z(1)=cos(x(2,1))*x(1,4)-x(1,1)*sin(x(2,1))*x(2,4) !d(x)/d(s2)
3102  z(2)=sin(x(2,1))*x(1,4)+x(1,1)*cos(x(2,1))*x(2,4) !d(y)/d(s2)
3103  CASE(3)
3104  z(1)=cos(x(2,1))*x(1,4)-x(1,1)*sin(x(2,1))*x(2,4) !d(x)/d(s2)
3105  z(2)=sin(x(2,1))*x(1,4)+x(1,1)*cos(x(2,1))*x(2,4) !d(y)/d(s2)
3106  z(3)=x(3,4) !d(z)/d(s2)
3107  CASE DEFAULT
3108  CALL flagerror("Invalid number of coordinates",err,error,*999)
3109  END SELECT
3110  ELSE
3111  CALL flagerror("Not enough X derivatives supplied",err,error,*999)
3112  ENDIF
3113  CASE(part_deriv_s2_s2)
3114  CALL flagerror("Not implemented",err,error,*999)
3115  CASE(part_deriv_s1_s2)
3116  IF(SIZE(x,2)>=6) THEN
3117  SELECT CASE(coordinate_system%NUMBER_OF_DIMENSIONS)
3118  CASE(2)
3119  z(1)=cos(x(2,1))*x(1,6)-x(1,2)*sin(x(2,1))*x(2,4)-&
3120  & (x(1,4)*sin(x(2,1))*x(2,2)+x(1,1)*cos(x(2,1))*&
3121  & x(2,2)*x(2,4)+x(1,1)*sin(x(2,1))*x(2,6)) !d2(x)/d(s1)d(s2)
3122  z(2)=sin(x(2,1))*x(1,6)+x(1,2)*cos(x(2,1))*x(2,4)+&
3123  & (x(1,4)*cos(x(2,1))*x(2,2)-x(1,1)*sin(x(2,1))*&
3124  & x(2,2)*x(2,4)+x(1,1)*cos(x(2,1))*x(2,6)) !d2(y)/d(s1)d(s2)
3125  CASE(3)
3126  z(1)=cos(x(2,1))*x(1,6)-x(1,2)*sin(x(2,1))*x(2,4)-&
3127  & (x(1,4)*sin(x(2,1))*x(2,2)+x(1,1)*cos(x(2,1))*&
3128  & x(2,2)*x(2,4)+x(1,1)*sin(x(2,1))*x(2,6)) !d2(x)/d(s1)d(s2)
3129  z(2)=sin(x(2,1))*x(1,6)+x(1,2)*cos(x(2,1))*x(2,4)+&
3130  & (x(1,4)*cos(x(2,1))*x(2,2)-x(1,1)*sin(x(2,1))*&
3131  & x(2,2)*x(2,4)+x(1,1)*cos(x(2,1))*x(2,6)) !d2(y)/d(s1)d(s2)
3132  z(3)=x(3,6) !d2(z)/d(s1)d(s2)
3133  CASE DEFAULT
3134  CALL flagerror("Invalid number of coordinates",err,error,*999)
3135  END SELECT
3136  ELSE
3137  CALL flagerror("Not enough X derivatives supplied",err,error,*999)
3138  ENDIF
3139  CASE(part_deriv_s3)
3140  IF(SIZE(x,2)>=7) THEN
3141  SELECT CASE(coordinate_system%NUMBER_OF_DIMENSIONS)
3142  CASE(2)
3143  z(1)=0.0_sp
3144  z(2)=0.0_sp
3145  CASE(3)
3146  z(1)=cos(x(2,1))*x(1,7)-x(1,1)*sin(x(2,1))*x(2,7) !d(x)/d(s3)
3147  z(2)=sin(x(2,1))*x(1,7)+x(1,1)*cos(x(2,1))*x(2,7) !d(y)/d(s3)
3148  z(3)=x(3,7) !d(z)/d(s3)
3149  CASE DEFAULT
3150  CALL flagerror("Invalid number of coordinates",err,error,*999)
3151  END SELECT
3152  ELSE
3153  CALL flagerror("Not enough X derivatives supplied",err,error,*999)
3154  ENDIF
3155  CASE(part_deriv_s3_s3)
3156  CALL flagerror("Not implemented",err,error,*999)
3157  CASE(part_deriv_s1_s3)
3158  IF(SIZE(x,2)>=9) THEN
3159  SELECT CASE(coordinate_system%NUMBER_OF_DIMENSIONS)
3160  CASE(2)
3161  z(1)=0.0_sp
3162  z(2)=0.0_sp
3163  CASE(3)
3164  z(1)=cos(x(2,1))*x(1,9)-x(1,2)*sin(x(2,1))*x(2,7)-&
3165  & (x(1,7)*sin(x(2,1))*x(2,2)+x(1,1)*cos(x(2,1))*&
3166  & x(2,2)*x(2,7)+x(1,1)*sin(x(2,1))*x(2,9)) !d2(x)/d(s1)d(s3)
3167  z(2)=sin(x(2,1))*x(1,9)+x(1,2)*cos(x(2,1))*x(2,7)+&
3168  & (x(1,7)*cos(x(2,1))*x(2,2)-x(1,1)*sin(x(2,1))*&
3169  & x(2,2)*x(2,7)+x(1,1)*cos(x(2,1))*x(2,9)) !d2(y)/d(s1)d(s3)
3170  z(3)=x(3,9) !d2(z)/d(s1)d(s3)
3171  CASE DEFAULT
3172  CALL flagerror("Invalid number of coordinates",err,error,*999)
3173  END SELECT
3174  ELSE
3175  CALL flagerror("Not enough X derivatives supplied",err,error,*999)
3176  ENDIF
3177  CASE(part_deriv_s2_s3)
3178  IF(SIZE(x,2)>=10) THEN
3179  SELECT CASE(coordinate_system%NUMBER_OF_DIMENSIONS)
3180  CASE(2)
3181  z(1)=0.0_sp
3182  z(2)=0.0_sp
3183  CASE(3)
3184  z(1)=cos(x(2,1))*x(1,10)-x(1,4)*sin(x(2,1))*x(2,7)-&
3185  & (x(1,7)*sin(x(2,1))*x(2,4)+x(1,1)*cos(x(2,1))*&
3186  & x(2,4)*x(2,7)+x(1,1)*sin(x(2,1))*x(2,10)) !d2(x)/d(s2)d(s3)
3187  z(2)=sin(x(2,1))*x(1,10)+x(1,4)*cos(x(2,1))*x(2,7)+&
3188  & (x(1,7)*cos(x(2,1))*x(2,4)-x(1,1)*sin(x(2,1))*&
3189  & x(2,4)*x(2,7)+x(1,1)*cos(x(2,1))*x(2,10)) !d2(y)/d(s2)d(s3)
3190  z(3)=x(3,10) !d2(z)/d(s2)d(s3)
3191  CASE DEFAULT
3192  CALL flagerror("Invalid number of coordinates",err,error,*999)
3193  END SELECT
3194  ELSE
3195  CALL flagerror("Not enough X derivatives supplied",err,error,*999)
3196  ENDIF
3197  CASE(part_deriv_s1_s2_s3)
3198  IF(SIZE(x,2)>=11) THEN
3199  SELECT CASE(coordinate_system%NUMBER_OF_DIMENSIONS)
3200  CASE(2)
3201  z(1)=0.0_sp
3202  z(2)=0.0_sp
3203  CASE(3)
3204  z(1)=-cos(x(2,1))*x(2,2)*x(2,4)*x(1,7)-&
3205  & sin(x(2,1))*(x(2,6)*x(1,7)+x(2,4)*x(1,9))-&
3206  & sin(x(2,1))*x(2,2)*x(1,10)+cos(x(2,1))*x(1,11)-&
3207  & cos(x(2,1))*x(2,2)*x(1,4)*x(2,7)-&
3208  & sin(x(2,1))*(x(1,6)*x(2,7)+x(1,4)*x(2,9))+&
3209  & (-cos(x(2,1))*x(1,2)+sin(x(2,1))*x(1,1)*x(2,2))*&
3210  & x(2,4)*x(2,7)-x(1,1)*cos(x(2,1))*(x(2,6)*x(2,7)+&
3211  & x(2,4)*x(2,9))+(-sin(x(2,1))*x(1,2)-x(1,1)*&
3212  & cos(x(2,1))*x(2,2))*x(2,10)-x(1,1)*&
3213  & sin(x(2,1))*x(2,11) !d3(x)/d(s1)d(s2)d(s3)
3214  z(2)=-sin(x(2,1))*x(2,2)*x(2,4)*x(1,7)+&
3215  & cos(x(2,1))*(x(2,6)*x(1,7)+x(2,4)*x(1,9))+&
3216  & cos(x(2,1))*x(2,2)*x(1,10)+sin(x(2,1))*x(1,11)-&
3217  & sin(x(2,1))*x(2,2)*x(1,4)*x(2,7)+&
3218  & cos(x(2,1))*(x(1,6)*x(2,7)+x(1,4)*x(2,9))+&
3219  & (-sin(x(2,1))*x(1,2)-x(1,1)*cos(x(2,1))*x(2,2))*&
3220  & x(2,4)*x(2,7)-x(1,1)*sin(x(2,1))*(x(2,6)*x(2,7)+&
3221  & x(2,4)*x(2,9))+(cos(x(2,1))*x(1,2)-x(1,1)*&
3222  & sin(x(2,1))*x(2,2))*x(2,10)+x(1,1)*&
3223  & cos(x(2,1))*x(2,11) !d3(y)/d(s1)d(s2)d(s3)
3224  z(3)=x(3,11) !d3(z)/d(s1)d(s2)d(s3)
3225  CASE DEFAULT
3226  CALL flagerror("Invalid number of coordinates",err,error,*999)
3227  END SELECT
3228  ELSE
3229  CALL flagerror("Not enough X derivatives supplied",err,error,*999)
3230  ENDIF
3231  CASE DEFAULT
3232  CALL flagerror("Invalid partial derivative type",err,error,*999)
3233  END SELECT
3234  CASE(coordinate_spherical_polar_type)
3235  IF(coordinate_system%NUMBER_OF_DIMENSIONS==3) THEN
3236  SELECT CASE(part_deriv_type)
3237  CASE(no_part_deriv)
3238  z=coordinate_convert_to_rc(coordinate_system,x(:,1),err,error)
3239  IF(err/=0) GOTO 999
3240  CASE(part_deriv_s1,part_deriv_s1_s1,part_deriv_s2,part_deriv_s2_s2,&
3241  & part_deriv_s1_s2,part_deriv_s3,part_deriv_s3_s3,&
3242  & part_deriv_s1_s3,part_deriv_s2_s3,part_deriv_s1_s2_s3)
3243  CALL flagerror("Not implemented",err,error,*999)
3244  CASE DEFAULT
3245  CALL flagerror("Invalid partial derivative type",err,error,*999)
3246  END SELECT
3247  ELSE
3248  CALL flagerror("Invalid number of coordinates",err,error,*999)
3249  ENDIF
3250  CASE(coordinate_prolate_spheroidal_type)
3251  IF(coordinate_system%NUMBER_OF_DIMENSIONS==3) THEN
3252  focus=REAL(coordinate_system%focus,sp)
3253  SELECT CASE(part_deriv_type)
3254  CASE(no_part_deriv)
3255  z=coordinate_convert_to_rc(coordinate_system,x(:,1),err,error)
3256  IF(err/=0) GOTO 999
3257  CASE(part_deriv_s1)
3258  IF(SIZE(x,2)>=2) THEN
3259  z(1)=focus*sinh(x(1,1))*cos(x(2,1))*x(1,2)-&
3260  & focus*cosh(x(1,1))*sin(x(2,1))*x(2,2) !d(x)/d(s1)
3261  z(2)=focus*cosh(x(1,1))*sin(x(2,1))*cos(x(3,1))*x(1,2)+&
3262  & focus*sinh(x(1,1))*cos(x(2,1))*cos(x(3,1))*x(2,2)-&
3263  & focus*sinh(x(1,1))*sin(x(2,1))*sin(x(3,1))*&
3264  & x(3,2) !d(y)/d(s1)
3265  z(3)=focus*cosh(x(1,1))*sin(x(2,1))*sin(x(3,1))*x(1,2)+&
3266  & focus*sinh(x(1,1))*cos(x(2,1))*sin(x(3,1))*x(2,2)+&
3267  & focus*sinh(x(1,1))*sin(x(2,1))*cos(x(3,1))*&
3268  & x(3,2) !d(z)/d(s1)
3269  ELSE
3270  CALL flagerror("Not enough X derivatives supplied", &
3271  & err,error,*999)
3272  ENDIF
3273  CASE(part_deriv_s1_s1)
3274  CALL flagerror("Not implemented",err,error,*999)
3275  CASE(part_deriv_s2)
3276  IF(SIZE(x,2)>=4) THEN
3277  z(1)=focus*sinh(x(1,1))*cos(x(2,1))*x(1,4)-&
3278  & focus*cosh(x(1,1))*sin(x(2,1))*x(2,4) !d(x)/d(s2)
3279  z(2)=focus*cosh(x(1,1))*sin(x(2,1))*cos(x(3,1))*x(1,4)+&
3280  & focus*sinh(x(1,1))*cos(x(2,1))*cos(x(3,1))*x(2,4)-&
3281  & focus*sinh(x(1,1))*sin(x(2,1))*sin(x(3,1))*&
3282  & x(3,4) !d(y)/d(s2)
3283  z(3)=focus*cosh(x(1,1))*sin(x(2,1))*sin(x(3,1))*x(1,4)+&
3284  & focus*sinh(x(1,1))*cos(x(2,1))*sin(x(3,1))*x(2,4)+&
3285  & focus*sinh(x(1,1))*sin(x(2,1))*cos(x(3,1))*&
3286  & x(3,4) !d(z)/d(s2)
3287  ELSE
3288  CALL flagerror("Not enough X derivatives supplied",&
3289  & err,error,*999)
3290  ENDIF
3291  CASE(part_deriv_s2_s2)
3292  CALL flagerror("Not implemented",err,error,*999)
3293  CASE(part_deriv_s1_s2)
3294  IF(SIZE(x,2)>=6) THEN
3295  z(1)=focus*(sinh(x(1,1))*cos(x(2,1))*x(1,6)+&
3296  & x(1,2)*(cosh(x(1,1))*cos(x(2,1))*x(1,4)-&
3297  & sinh(x(1,1))*sin(x(2,1))*x(2,4))-&
3298  & cosh(x(1,1))*sin(x(2,1))*x(2,6)-&
3299  & x(2,2)*(sinh(x(1,1))*sin(x(2,1))*x(1,4)+&
3300  & cosh(x(1,1))*cos(x(2,1))*x(2,4))) !d2(x)/d(s1)d(s2)
3301  z(2)=focus*(cosh(x(1,1))*sin(x(2,1))*cos(x(3,1))*x(1,6)+&
3302  & x(1,2)*(sinh(x(1,1))*sin(x(2,1))*cos(x(3,1))*&
3303  & x(1,4)+cosh(x(1,1))*cos(x(2,1))*cos(x(3,1))*&
3304  & x(2,4)-cosh(x(1,1))*sin(x(2,1))*sin(x(3,1))*&
3305  & x(3,4))+sinh(x(1,1))*cos(x(2,1))*cos(x(3,1))*&
3306  & x(2,6)+x(2,2)*(cosh(x(1,1))*cos(x(2,1))*cos(x(3,1))*&
3307  & x(1,4)-sinh(x(1,1))*sin(x(2,1))*cos(x(3,1))*x(2,4)-&
3308  & sinh(x(1,1))*cos(x(2,1))*sin(x(3,1))*x(3,4))-&
3309  & sinh(x(1,1))*sin(x(2,1))*sin(x(3,1))*x(3,6)-&
3310  & x(3,2)*(cosh(x(1,1))*sin(x(2,1))*sin(x(3,1))*&
3311  & x(1,4)+sinh(x(1,1))*cos(x(2,1))*sin(x(3,1))*&
3312  & x(2,4)+sinh(x(1,1))*sin(x(2,1))*cos(x(3,1))*&
3313  & x(3,4))) !d2(y)/d(s1)d(s2)
3314  z(3)=focus*(cosh(x(1,1))*sin(x(2,1))*sin(x(3,1))*x(1,6)+&
3315  & x(1,2)*(sinh(x(1,1))*sin(x(2,1))*sin(x(3,1))*&
3316  & x(1,4)+cosh(x(1,1))*cos(x(2,1))*sin(x(3,1))*&
3317  & x(2,4)+cosh(x(1,1))*sin(x(2,1))*cos(x(3,1))*&
3318  & x(3,4))+sinh(x(1,1))*cos(x(2,1))*sin(x(3,1))*x(2,6)+&
3319  & x(2,2)*(cosh(x(1,1))*cos(x(2,1))*sin(x(3,1))*&
3320  & x(1,4)-sinh(x(1,1))*sin(x(2,1))*sin(x(3,1))*&
3321  & x(2,4)+sinh(x(1,1))*cos(x(2,1))*cos(x(3,1))*&
3322  & x(3,4))+sinh(x(1,1))*sin(x(2,1))*cos(x(3,1))*x(3,6)+&
3323  & x(3,2)*(cosh(x(1,1))*sin(x(2,1))*cos(x(3,1))*&
3324  & x(1,4)+sinh(x(1,1))*cos(x(2,1))*cos(x(3,1))*&
3325  & x(2,4)-sinh(x(1,1))*sin(x(2,1))*sin(x(3,1))*&
3326  & x(3,4))) !d2(z)/d(s1)d(s2)
3327  ELSE
3328  CALL flagerror("Not enough X derivatives supplied",&
3329  & err,error,*999)
3330  ENDIF
3331  CASE(part_deriv_s3)
3332  IF(SIZE(x,2)>=7) THEN
3333  z(1)=focus*sinh(x(1,1))*cos(x(2,1))*x(1,7)-&
3334  & focus*cosh(x(1,1))*sin(x(2,1))*x(2,7) !d(x)/d(s3)
3335  z(2)=focus*cosh(x(1,1))*sin(x(2,1))*cos(x(3,1))*x(1,7)+&
3336  & focus*sinh(x(1,1))*cos(x(2,1))*cos(x(3,1))*x(2,7)-&
3337  & focus*sinh(x(1,1))*sin(x(2,1))*sin(x(3,1))*&
3338  & x(3,7) !d(y)/d(s3)
3339  z(3)=focus*cosh(x(1,1))*sin(x(2,1))*sin(x(3,1))*x(1,7)+&
3340  & focus*sinh(x(1,1))*cos(x(2,1))*sin(x(3,1))*x(2,7)+&
3341  & focus*sinh(x(1,1))*sin(x(2,1))*cos(x(3,1))*&
3342  & x(3,7) !d(z)/d(s3)
3343  ELSE
3344  CALL flagerror("Not enough X derivatives supplied",&
3345  & err,error,*999)
3346  ENDIF
3347  CASE(part_deriv_s3_s3)
3348  CALL flagerror("Not implemented",err,error,*999)
3349  CASE(part_deriv_s1_s3)
3350  IF(SIZE(x,2)>=9) THEN
3351  z(1)=focus*(sinh(x(1,1))*cos(x(2,1))*x(1,9)+&
3352  & x(1,2)*(cosh(x(1,1))*cos(x(2,1))*x(1,7)-&
3353  & sinh(x(1,1))*sin(x(2,1))*x(2,7))-&
3354  & cosh(x(1,1))*sin(x(2,1))*x(2,9)-&
3355  & x(2,2)*(sinh(x(1,1))*sin(x(2,1))*x(1,7)+&
3356  & cosh(x(1,1))*cos(x(2,1))*x(2,7))) !d2(x)/d(s1)d(s3)
3357  z(2)=focus*(cosh(x(1,1))*sin(x(2,1))*cos(x(3,1))*x(1,9)+&
3358  & x(1,2)*(sinh(x(1,1))*sin(x(2,1))*cos(x(3,1))*x(1,7)+&
3359  & cosh(x(1,1))*cos(x(2,1))*cos(x(3,1))*x(2,7)-&
3360  & cosh(x(1,1))*sin(x(2,1))*sin(x(3,1))*x(3,7))+&
3361  & sinh(x(1,1))*cos(x(2,1))*cos(x(3,1))*x(2,9)+&
3362  & x(2,2)*(cosh(x(1,1))*cos(x(2,1))*cos(x(3,1))*x(1,7)-&
3363  & sinh(x(1,1))*sin(x(2,1))*cos(x(3,1))*x(2,7)-&
3364  & sinh(x(1,1))*cos(x(2,1))*sin(x(3,1))*x(3,7))-&
3365  & sinh(x(1,1))*sin(x(2,1))*sin(x(3,1))*x(3,9)-&
3366  & x(3,2)*(cosh(x(1,1))*sin(x(2,1))*sin(x(3,1))*x(1,7)+&
3367  & sinh(x(1,1))*cos(x(2,1))*sin(x(3,1))* x(2,7)+&
3368  & sinh(x(1,1))*sin(x(2,1))*cos(x(3,1))*&
3369  & x(3,7))) !d2(y)/d(s1)d(s3)
3370  z(3)=focus*(cosh(x(1,1))*sin(x(2,1))*sin(x(3,1))*x(1,9)+&
3371  & x(1,2)*(sinh(x(1,1))*sin(x(2,1))*sin(x(3,1))*x(1,7)+&
3372  & cosh(x(1,1))*cos(x(2,1))*sin(x(3,1))*x(2,7)+&
3373  & cosh(x(1,1))*sin(x(2,1))*cos(x(3,1))*x(3,7))+&
3374  & sinh(x(1,1))*cos(x(2,1))*sin(x(3,1))*x(2,9)+&
3375  & x(2,2)*(cosh(x(1,1))*cos(x(2,1))*sin(x(3,1))*x(1,7)-&
3376  & sinh(x(1,1))*sin(x(2,1))*sin(x(3,1))*x(2,7)+&
3377  & sinh(x(1,1))*cos(x(2,1))*cos(x(3,1))*x(3,7))+&
3378  & sinh(x(1,1))*sin(x(2,1))*cos(x(3,1))*x(3,9)+&
3379  & x(3,2)*(cosh(x(1,1))*sin(x(2,1))*cos(x(3,1))*x(1,7)+&
3380  & sinh(x(1,1))*cos(x(2,1))*cos(x(3,1))*x(2,7)-&
3381  & sinh(x(1,1))*sin(x(2,1))*sin(x(3,1))*&
3382  & x(3,7))) !d2(z)/d(s1)d(s3)
3383  ELSE
3384  CALL flagerror("Not enough X derivatives supplied",&
3385  & err,error,*999)
3386  ENDIF
3387  CASE(part_deriv_s2_s3)
3388  IF(SIZE(x,2)>=10) THEN
3389  z(1)=focus*(sinh(x(1,1))*cos(x(2,1))*x(1,10)+&
3390  & x(1,4)*(cosh(x(1,1))*cos(x(2,1))*x(1,7)-&
3391  & sinh(x(1,1))*sin(x(2,1))*x(2,7))-&
3392  & cosh(x(1,1))*sin(x(2,1))*x(2,10)-&
3393  & x(2,4)*(sinh(x(1,1))*sin(x(2,1))*x(1,7)+&
3394  & cosh(x(1,1))*cos(x(2,1))*x(2,7))) !d2(x)/d(s2)d(s3)
3395  z(2)=focus*(cosh(x(1,1))*sin(x(2,1))*cos(x(3,1))*x(1,10)+&
3396  & x(1,4)*(sinh(x(1,1))*sin(x(2,1))*cos(x(3,1))*x(1,7)+&
3397  & cosh(x(1,1))*cos(x(2,1))*cos(x(3,1))*x(2,7)-&
3398  & cosh(x(1,1))*sin(x(2,1))*sin(x(3,1))*x(3,7))+&
3399  & sinh(x(1,1))*cos(x(2,1))*cos(x(3,1))*x(2,10)+&
3400  & x(2,4)*(cosh(x(1,1))*cos(x(2,1))*cos(x(3,1))*x(1,7)-&
3401  & sinh(x(1,1))*sin(x(2,1))*cos(x(3,1))*x(2,7)-&
3402  & sinh(x(1,1))*cos(x(2,1))*sin(x(3,1))*x(3,7))-&
3403  & sinh(x(1,1))*sin(x(2,1))*sin(x(3,1))*x(3,10)-&
3404  & x(3,4)*(cosh(x(1,1))*sin(x(2,1))*sin(x(3,1))*x(1,7)+&
3405  & sinh(x(1,1))*cos(x(2,1))*sin(x(3,1))*x(2,7)+&
3406  & sinh(x(1,1))*sin(x(2,1))*cos(x(3,1))*&
3407  & x(3,7))) !d2(y)/d(s2)d(s3)
3408  z(3)=focus*(cosh(x(1,1))*sin(x(2,1))*sin(x(3,1))*x(1,10)+&
3409  & x(1,4)*(sinh(x(1,1))*sin(x(2,1))*sin(x(3,1))*x(1,7)+&
3410  & cosh(x(1,1))*cos(x(2,1))*sin(x(3,1))*x(2,7)+&
3411  & cosh(x(1,1))*sin(x(2,1))*cos(x(3,1))*x(3,7))+&
3412  & sinh(x(1,1))*cos(x(2,1))*sin(x(3,1))*x(2,10)+&
3413  & x(2,4)*(cosh(x(1,1))*cos(x(2,1))*sin(x(3,1))*x(1,7)-&
3414  & sinh(x(1,1))*sin(x(2,1))*sin(x(3,1))*x(2,7)+&
3415  & sinh(x(1,1))*cos(x(2,1))*cos(x(3,1))*x(3,7))+&
3416  & sinh(x(1,1))*sin(x(2,1))*cos(x(3,1))*x(3,10)+&
3417  & x(3,4)*(cosh(x(1,1))*sin(x(2,1))*cos(x(3,1))*x(1,7)+&
3418  & sinh(x(1,1))*cos(x(2,1))*cos(x(3,1))*x(2,7)-&
3419  & sinh(x(1,1))*sin(x(2,1))*sin(x(3,1))*&
3420  & x(3,7))) !d2(z)/d(s2)d(s3)
3421  ELSE
3422  CALL flagerror("Not enough X derivatives supplied",&
3423  & err,error,*999)
3424  ENDIF
3425  CASE(part_deriv_s1_s2_s3)
3426  IF(SIZE(x,2)>=11) THEN
3427  z(1)=focus*((sinh(x(1,1))*cos(x(2,1))*x(1,2)-&
3428  & cosh(x(1,1))*sin(x(2,1))*x(2,2))*x(1,4)*x(1,7)+&
3429  & cosh(x(1,1))*cos(x(2,1))*(x(1,6)*x(1,7)+x(1,4)*x(1,9))+&
3430  & (-cosh(x(1,1))*sin(x(2,1))*x(1,2)-&
3431  & sinh(x(1,1))*cos(x(2,1))*x(2,2))*x(2,4)*x(1,7)-&
3432  & sinh(x(1,1))*sin(x(2,1))*(x(2,6)*x(1,7)+x(2,4)*x(1,9))+&
3433  & (cosh(x(1,1))*cos(x(2,1))*x(1,2)-&
3434  & sinh(x(1,1))*sin(x(2,1))*x(2,2))*x(1,10)+&
3435  & sinh(x(1,1))*cos(x(2,1))*x(1,11)+&
3436  & (-cosh(x(1,1))*sin(x(2,1))*x(1,2)-&
3437  & sinh(x(1,1))*cos(x(2,1))*x(2,2))*x(1,4)*x(2,7)-&
3438  & sinh(x(1,1))*sin(x(2,1))*(x(1,6)*x(2,7)+x(1,4)*x(2,9))+&
3439  & (-sinh(x(1,1))*cos(x(2,1))*x(1,2)+&
3440  & cosh(x(1,1))*sin(x(2,1))*x(2,2))*x(2,4)*x(2,7)-&
3441  & cosh(x(1,1))*cos(x(2,1))*(x(2,6)*x(2,7)+x(2,4)*x(2,9))+&
3442  & (-sinh(x(1,1))*sin(x(2,1))*x(1,2)-&
3443  & cosh(x(1,1))*cos(x(2,1))*x(2,2))*x(2,10)-&
3444  & cosh(x(1,1))*sin(x(2,1))*x(2,11)) !d3(x)/d(s1)d(s2)d(s3)
3445  z(2)=focus*((cosh(x(1,1))*sin(x(2,1))*cos(x(3,1))*x(1,2)+&
3446  & sinh(x(1,1))*cos(x(2,1))*cos(x(3,1))*x(2,2)-&
3447  & sinh(x(1,1))*sin(x(2,1))*sin(x(3,1))*x(3,2))*&
3448  & x(1,4)*x(1,7)+sinh(x(1,1))*sin(x(2,1))*cos(x(3,1))*&
3449  & (x(1,6)*x(1,7)+x(1,4)*x(1,9))+&
3450  & (sinh(x(1,1))*cos(x(2,1))*cos(x(3,1))*x(1,2)-&
3451  & cosh(x(1,1))*sin(x(2,1))*cos(x(3,1))*x(2,2)-&
3452  & cosh(x(1,1))*cos(x(2,1))*sin(x(3,1))*x(3,2))*&
3453  & x(2,4)*x(1,7)+cosh(x(1,1))*cos(x(2,1))*cos(x(3,1))*&
3454  & (x(2,6)*x(1,7)+x(2,4)*x(1,9))+&
3455  & (-sinh(x(1,1))*sin(x(2,1))*sin(x(3,1))*x(1,2)-&
3456  & cosh(x(1,1))*cos(x(2,1))*sin(x(3,1))*x(2,2)-&
3457  & cosh(x(1,1))*sin(x(2,1))*cos(x(3,1))*x(3,2))*&
3458  & x(3,4)*x(1,7)-cosh(x(1,1))*sin(x(2,1))*sin(x(3,1))*&
3459  & (x(3,6)*x(1,7)+x(3,4)*x(1,9))+&
3460  & (sinh(x(1,1))*sin(x(2,1))*cos(x(3,1))*x(1,2)+&
3461  & cosh(x(1,1))*cos(x(2,1))*cos(x(3,1))*x(2,2)-&
3462  & cosh(x(1,1))*sin(x(2,1))*sin(x(3,1))*x(3,2))*x(1,10)+&
3463  & cosh(x(1,1))*sin(x(2,1))*cos(x(3,1))*x(1,11))+&
3464  & focus*((sinh(x(1,1))*cos(x(2,1))*cos(x(3,1))*x(1,2)-&
3465  & cosh(x(1,1))*sin(x(2,1))*cos(x(3,1))*x(2,2)-&
3466  & cosh(x(1,1))*cos(x(2,1))*sin(x(3,1))*x(3,2))*&
3467  & x(1,4)*x(2,7)+cosh(x(1,1))*cos(x(2,1))*cos(x(3,1))*&
3468  & (x(1,6)*x(2,7)+x(1,4)*x(2,9))+&
3469  & (-cosh(x(1,1))*sin(x(2,1))*cos(x(3,1))*x(1,2)-&
3470  & sinh(x(1,1))*cos(x(2,1))*cos(x(3,1))*x(2,2)+&
3471  & sinh(x(1,1))*sin(x(2,1))*sin(x(3,1))*x(3,2))*&
3472  & x(2,4)*x(2,7)-sinh(x(1,1))*sin(x(2,1))*cos(x(3,1))*&
3473  & (x(2,6)*x(2,7)+x(2,4)*x(2,9))+&
3474  & (-cosh(x(1,1))*cos(x(2,1))*sin(x(3,1))*x(1,2)+&
3475  & sinh(x(1,1))*sin(x(2,1))*sin(x(3,1))*x(2,2)-&
3476  & sinh(x(1,1))*cos(x(2,1))*cos(x(3,1))*x(3,2))*&
3477  & x(3,4)*x(2,7)-sinh(x(1,1))*cos(x(2,1))*sin(x(3,1))*&
3478  & (x(3,6)*x(2,7)+x(3,4)*x(2,9))+&
3479  & (cosh(x(1,1))*cos(x(2,1))*cos(x(3,1))*x(1,2)-&
3480  & sinh(x(1,1))*sin(x(2,1))*cos(x(3,1))*x(2,2)-&
3481  & sinh(x(1,1))*cos(x(2,1))*sin(x(3,1))*x(3,2))*x(2,10)+&
3482  & sinh(x(1,1))*cos(x(2,1))*cos(x(3,1))*x(2,11))+&
3483  & focus*((-sinh(x(1,1))*sin(x(2,1))*sin(x(3,1))*x(1,2)-&
3484  & cosh(x(1,1))*cos(x(2,1))*sin(x(3,1))*x(2,2)-&
3485  & cosh(x(1,1))*sin(x(2,1))*cos(x(3,1))*x(3,2))*&
3486  & x(1,4)*x(3,7)-cosh(x(1,1))*sin(x(2,1))*sin(x(3,1))*&
3487  & (x(1,6)*x(3,7)+x(1,4)*x(3,9))+&
3488  & (-cosh(x(1,1))*cos(x(2,1))*sin(x(3,1))*x(1,2)+&
3489  & sinh(x(1,1))*sin(x(2,1))*sin(x(3,1))*x(2,2)-&
3490  & sinh(x(1,1))*cos(x(2,1))*cos(x(3,1))*x(3,2))*&
3491  & x(2,4)*x(3,7)-sinh(x(1,1))*cos(x(2,1))*sin(x(3,1))*&
3492  & (x(2,6)*x(3,7)+x(2,4)*x(3,9))+&
3493  & (-cosh(x(1,1))*sin(x(2,1))*cos(x(3,1))*x(1,2)-&
3494  & sinh(x(1,1))*cos(x(2,1))*cos(x(3,1))*x(2,2)+&
3495  & sinh(x(1,1))*sin(x(2,1))*sin(x(3,1))*x(3,2))*&
3496  & x(3,4)*x(3,7)-sinh(x(1,1))*sin(x(2,1))*cos(x(3,1))*&
3497  & (x(3,6)*x(3,7)+x(3,4)*x(3,9))+&
3498  & (-cosh(x(1,1))*sin(x(2,1))*sin(x(3,1))*x(1,2)-&
3499  & sinh(x(1,1))*cos(x(2,1))*sin(x(3,1))*x(2,2)-&
3500  & sinh(x(1,1))*sin(x(2,1))*cos(x(3,1))*x(3,2))*x(3,10)-&
3501  & sinh(x(1,1))*sin(x(2,1))*sin(x(3,1))*&
3502  & x(3,11)) !d3(y)/d(s1)d(s2)d(s3)
3503  z(3)=focus*((cosh(x(1,1))*sin(x(2,1))*sin(x(3,1))*x(1,2)+&
3504  & sinh(x(1,1))*cos(x(2,1))*sin(x(3,1))*x(2,2)+&
3505  & sinh(x(1,1))*sin(x(2,1))*cos(x(3,1))*x(3,2))*&
3506  & x(1,4)*x(1,7)+sinh(x(1,1))*sin(x(2,1))*sin(x(3,1))*&
3507  & (x(1,6)*x(1,7)+x(1,4)*x(1,9))+&
3508  & (sinh(x(1,1))*cos(x(2,1))*sin(x(3,1))*x(1,2)-&
3509  & cosh(x(1,1))*sin(x(2,1))*sin(x(3,1))*x(2,2)+&
3510  & cosh(x(1,1))*cos(x(2,1))*cos(x(3,1))*x(3,2))*&
3511  & x(2,4)*x(1,7)+cosh(x(1,1))*cos(x(2,1))*sin(x(3,1))*&
3512  & (x(2,6)*x(1,7)+x(2,4)*x(1,9))+&
3513  & (sinh(x(1,1))*sin(x(2,1))*cos(x(3,1))*x(1,2)+&
3514  & cosh(x(1,1))*cos(x(2,1))*cos(x(3,1))*x(2,2)-&
3515  & cosh(x(1,1))*sin(x(2,1))*sin(x(3,1))*x(3,2))*&
3516  & x(3,4)*x(1,7)+cosh(x(1,1))*sin(x(2,1))*cos(x(3,1))*&
3517  & (x(3,6)*x(1,7)+x(3,4)*x(1,9))+&
3518  & (sinh(x(1,1))*sin(x(2,1))*sin(x(3,1))*x(1,2)+&
3519  & cosh(x(1,1))*cos(x(2,1))*sin(x(3,1))*x(2,2)+&
3520  & cosh(x(1,1))*sin(x(2,1))*cos(x(3,1))*x(3,2))*x(1,10)+&
3521  & cosh(x(1,1))*sin(x(2,1))*sin(x(3,1))*x(1,11))+&
3522  & focus*((sinh(x(1,1))*cos(x(2,1))*sin(x(3,1))*x(1,2)-&
3523  & cosh(x(1,1))*sin(x(2,1))*sin(x(3,1))*x(2,2)+&
3524  & cosh(x(1,1))*cos(x(2,1))*cos(x(3,1))*x(3,2))*&
3525  & x(1,4)*x(2,7)+cosh(x(1,1))*cos(x(2,1))*sin(x(3,1))*&
3526  & (x(1,6)*x(2,7)+x(1,4)*x(2,9))+&
3527  & (-cosh(x(1,1))*sin(x(2,1))*sin(x(3,1))*x(1,2)-&
3528  & sinh(x(1,1))*cos(x(2,1))*sin(x(3,1))*x(2,2)-&
3529  & sinh(x(1,1))*sin(x(2,1))*cos(x(3,1))*x(3,2))*&
3530  & x(2,4)*x(2,7)-sinh(x(1,1))*sin(x(2,1))*sin(x(3,1))*&
3531  & (x(2,6)*x(2,7)+x(2,4)*x(2,9))+&
3532  & (cosh(x(1,1))*cos(x(2,1))*cos(x(3,1))*x(1,2)-&
3533  & sinh(x(1,1))*sin(x(2,1))*cos(x(3,1))*x(2,2)-&
3534  & sinh(x(1,1))*cos(x(2,1))*sin(x(3,1))*x(3,2))*&
3535  & x(3,4)*x(2,7)+sinh(x(1,1))*cos(x(2,1))*cos(x(3,1))*&
3536  & (x(3,6)*x(2,7)+x(3,4)*x(2,9))+&
3537  & (cosh(x(1,1))*cos(x(2,1))*sin(x(3,1))*x(1,2)-&
3538  & sinh(x(1,1))*sin(x(2,1))*sin(x(3,1))*x(2,2)+&
3539  & sinh(x(1,1))*cos(x(2,1))*cos(x(3,1))*x(3,2))*x(2,10)+&
3540  & sinh(x(1,1))*cos(x(2,1))*sin(x(3,1))*x(2,11))+&
3541  & focus*((sinh(x(1,1))*sin(x(2,1))*cos(x(3,1))*x(1,2)+&
3542  & cosh(x(1,1))*cos(x(2,1))*cos(x(3,1))*x(2,2)-&
3543  & cosh(x(1,1))*sin(x(2,1))*sin(x(3,1))*x(3,2))*&
3544  & x(1,4)*x(3,7)+cosh(x(1,1))*sin(x(2,1))*cos(x(3,1))*&
3545  & (x(1,6)*x(3,7)+x(1,4)*x(3,9))+&
3546  & (cosh(x(1,1))*cos(x(2,1))*cos(x(3,1))*x(1,2)-&
3547  & sinh(x(1,1))*sin(x(2,1))*cos(x(3,1))*x(2,2)-&
3548  & sinh(x(1,1))*cos(x(2,1))*sin(x(3,1))*x(3,2))*&
3549  & x(2,4)*x(3,7)+sinh(x(1,1))*cos(x(2,1))*cos(x(3,1))*&
3550  & (x(2,6)*x(3,7)+x(2,4)*x(3,9))+&
3551  & (-cosh(x(1,1))*sin(x(2,1))*sin(x(3,1))*x(1,2)-&
3552  & sinh(x(1,1))*cos(x(2,1))*sin(x(3,1))*x(2,2)-&
3553  & sinh(x(1,1))*sin(x(2,1))*cos(x(3,1))*x(3,2))*&
3554  & x(3,4)*x(3,7)-sinh(x(1,1))*sin(x(2,1))*sin(x(3,1))*&
3555  & (x(3,6)*x(3,7)+x(3,4)*x(3,9))+&
3556  & (cosh(x(1,1))*sin(x(2,1))*cos(x(3,1))*x(1,2)+&
3557  & sinh(x(1,1))*cos(x(2,1))*cos(x(3,1))*x(2,2)-&
3558  & sinh(x(1,1))*sin(x(2,1))*sin(x(3,1))*x(3,2))*x(3,10)+&
3559  & sinh(x(1,1))*sin(x(2,1))*cos(x(3,1))*&
3560  & x(3,11)) !d3(z)/d(s1)d(s2)d(s3)
3561  ELSE
3562  CALL flagerror("Not enough X derivatives supplied",&
3563  & err,error,*999)
3564  ENDIF
3565  CASE DEFAULT
3566  CALL flagerror("Invalid partial derivative type",err,error,*999)
3567  END SELECT
3568  ENDIF
3569  CASE(coordinate_oblate_spheroidal_type)
3570  IF(coordinate_system%NUMBER_OF_DIMENSIONS==3) THEN
3571  SELECT CASE(part_deriv_type)
3572  CASE(no_part_deriv)
3573  z=coordinate_convert_to_rc(coordinate_system,x(:,1),err,error)
3574  IF(err/=0) GOTO 999
3575  CASE(part_deriv_s1,part_deriv_s1_s1,part_deriv_s2,part_deriv_s2_s2,&
3576  & part_deriv_s1_s2,part_deriv_s3,part_deriv_s3_s3,&
3577  & part_deriv_s1_s3,part_deriv_s2_s3,part_deriv_s1_s2_s3)
3578  CALL flagerror("Not implemented",err,error,*999)
3579  CASE DEFAULT
3580  CALL flagerror("Invalid partial derivative type",err,error,*999)
3581  END SELECT
3582  ELSE
3583  CALL flagerror("Invalid number of coordinates",err,error,*999)
3584  ENDIF
3585  CASE DEFAULT
3586  CALL flagerror("Invalid coordinate type",err,error,*999)
3587  END SELECT
3588  ELSE
3589  CALL flagerror("The sizes of the vectors X and Z do not match",&
3590  & err,error,*999)
3591  ENDIF
3592 
3593  exits("COORDINATE_DERIVATIVE_CONVERT_TO_RC_SP")
3594  RETURN
3595 999 errorsexits("COORDINATE_DERIVATIVE_CONVERT_TO_RC_SP",err,error)
3596  RETURN 1
3597  END SUBROUTINE coordinate_derivative_convert_to_rc_sp
3598 
3599  !
3600  !================================================================================================================================
3601  !
3602 
3606  SUBROUTINE coordinate_derivative_norm(COORDINATE_SYSTEM,PART_DERIV_INDEX,INTERPOLATED_POINT,DERIV_NORM,ERR,ERROR,*)
3607 
3608  !Argument variables
3609  TYPE(coordinate_system_type), POINTER :: COORDINATE_SYSTEM
3610  INTEGER(INTG), INTENT(IN) :: PART_DERIV_INDEX
3611  TYPE(field_interpolated_point_type), POINTER :: INTERPOLATED_POINT
3612  REAL(DP), INTENT(OUT) :: DERIV_NORM
3613  INTEGER(INTG), INTENT(OUT) :: ERR
3614  TYPE(varying_string), INTENT(OUT) :: ERROR
3615  !Local variables
3616  INTEGER(INTG) :: component_idx,NUMBER_OF_COMPONENTS
3617  REAL(DP) :: FOCUS,SL,SM
3618  TYPE(varying_string) :: LOCAL_ERROR
3619 
3620  enters("COORDINATE_DERIVATIVE_NORM",err,error,*999)
3621 
3622  deriv_norm=0.0_dp
3623  IF(ASSOCIATED(coordinate_system)) THEN
3624  IF(ASSOCIATED(interpolated_point)) THEN
3625  IF(interpolated_point%PARTIAL_DERIVATIVE_TYPE>=first_part_deriv) THEN
3626  IF(part_deriv_index<=interpolated_point%MAX_PARTIAL_DERIVATIVE_INDEX) THEN
3627  number_of_components=SIZE(interpolated_point%VALUES,1)
3628  SELECT CASE(part_deriv_index)
3629  CASE(part_deriv_s1,part_deriv_s2,part_deriv_s3)
3630  SELECT CASE(coordinate_system%TYPE)
3631  CASE(coordinate_rectangular_cartesian_type)
3632  DO component_idx=1,number_of_components
3633  deriv_norm=deriv_norm+interpolated_point%VALUES(component_idx,part_deriv_index)**2
3634  ENDDO !component_index
3635  CASE(coordinate_cylindrical_polar_type)
3636  IF(number_of_components==2) THEN
3637  deriv_norm=interpolated_point%VALUES(1,part_deriv_index)**2+(interpolated_point%VALUES(1,1)* &
3638  & interpolated_point%VALUES(2,part_deriv_index))**2
3639  ELSE IF(number_of_components==3) THEN
3640  deriv_norm=interpolated_point%VALUES(1,part_deriv_index)**2+(interpolated_point%VALUES(1,1)* &
3641  & interpolated_point%VALUES(2,part_deriv_index))**2+interpolated_point%VALUES(3,part_deriv_index)**2
3642  ELSE
3643  local_error="The number of components for the interpolated point of "// &
3644  & trim(number_to_vstring(number_of_components,"*",err,error))// &
3645  & " is invalid for a cylindrical polar coordinate system."
3646  CALL flagerror(local_error,err,error,*999)
3647  ENDIF
3648  CASE(coordinate_spherical_polar_type)
3649  deriv_norm=interpolated_point%VALUES(1,part_deriv_index)**2+(interpolated_point%VALUES(1,1)* &
3650  & interpolated_point%VALUES(2,part_deriv_index)*cos(interpolated_point%VALUES(3,1)))**2+ &
3651  & (interpolated_point%VALUES(1,1)*interpolated_point%VALUES(3,part_deriv_index))**2
3652  CASE(coordinate_prolate_spheroidal_type)
3653  focus=coordinate_system%FOCUS
3654  sl=sinh(interpolated_point%VALUES(1,1))
3655  sm=sin(interpolated_point%VALUES(2,1))
3656  deriv_norm=focus*focus*((sl*sl+sm*sm)*(interpolated_point%VALUES(1,part_deriv_index)**2+ &
3657  & interpolated_point%VALUES(2,part_deriv_index))**2)+(sl*sm*interpolated_point%VALUES(3,part_deriv_index))**2
3658  CASE(coordinate_oblate_spheroidal_type)
3659  CALL flagerror("Not implemented.",err,error,*999)
3660  CASE DEFAULT
3661  local_error="The coordinate system type of "//trim(number_to_vstring(coordinate_system%TYPE,"*",err,error))// &
3662  & " is invalid."
3663  CALL flagerror(local_error,err,error,*999)
3664  END SELECT
3665  deriv_norm=sqrt(deriv_norm)
3666  CASE DEFAULT
3667  local_error="The partial derivative index of "//trim(number_to_vstring(part_deriv_index,"*",err,error))// &
3668  & " is invalid."
3669  CALL flagerror(local_error,err,error,*999)
3670  END SELECT
3671  ELSE
3672  local_error="The partial derivative index of "//trim(number_to_vstring(part_deriv_index,"*",err,error))// &
3673  & " is invalid. The interpolated point has a maximum number of partial derivatives of "// &
3674  & trim(number_to_vstring(interpolated_point%MAX_PARTIAL_DERIVATIVE_INDEX,"*",err,error))//"."
3675  CALL flagerror(local_error,err,error,*999)
3676  ENDIF
3677  ELSE
3678  CALL flagerror("The point has not been interpolated to include first derivative values.",err,error,*999)
3679  ENDIF
3680  ELSE
3681  CALL flagerror("Interpolated point is not associated.",err,error,*999)
3682  ENDIF
3683  ELSE
3684  CALL flagerror("Coordinate system is not associated.",err,error,*999)
3685  ENDIF
3686 
3687  exits("COORDINATE_DERIVATIVE_NORM")
3688  RETURN
3689 999 errorsexits("COORDINATE_DERIVATIVE_NORM",err,error)
3690  RETURN 1
3691  END SUBROUTINE coordinate_derivative_norm
3692 
3693  !
3694  !================================================================================================================================
3695  !
3696 
3698  SUBROUTINE coordinate_interpolation_adjust(COORDINATE_SYSTEM,PARTIAL_DERIVATIVE_INDEX,VALUE,ERR,ERROR,*)
3699 
3700  !Argument variables
3701  TYPE(coordinate_system_type), POINTER :: COORDINATE_SYSTEM
3702  INTEGER(INTG), INTENT(IN) :: PARTIAL_DERIVATIVE_INDEX
3703  REAL(DP), INTENT(INOUT) :: VALUE
3704  INTEGER(INTG), INTENT(OUT) :: ERR
3705  TYPE(varying_string), INTENT(OUT) :: ERROR
3706  !Local variables
3707  REAL(DP) :: COSHX,CSS,D,DES,FOCUS,R,SS,SINHX,THETA
3708  TYPE(varying_string) :: LOCAL_ERROR
3709 
3710  enters("COORDINATE_INTERPOLATION_ADJUST",err,error,*999)
3711 
3712  IF(ASSOCIATED(coordinate_system)) THEN
3713  SELECT CASE(coordinate_system%TYPE)
3714  CASE(coordinate_rectangular_cartesian_type)
3715  !Do nothing
3716  CASE(coordinate_cylindrical_polar_type)
3717  SELECT CASE(coordinate_system%RADIAL_INTERPOLATION_TYPE)
3718  CASE(coordinate_radial_interpolation_type)
3719  !Do nothing
3720  CASE(coordinate_radial_squared_interpolation_type)
3721  r=sqrt(VALUE)
3722  IF(partial_derivative_index==no_part_deriv) THEN
3723  VALUE=r
3724  ELSE
3725  VALUE=VALUE/(2.0_dp*r)
3726  ENDIF
3727  CASE DEFAULT
3728  local_error="The radial interpolation type of "//trim(number_to_vstring(coordinate_system% &
3729  & radial_interpolation_type,"*",err,error))//" is invalid for a cylindrical coordinate system."
3730  CALL flagerror(local_error,err,error,*999)
3731  END SELECT
3732  CASE(coordinate_spherical_polar_type)
3733  SELECT CASE(coordinate_system%RADIAL_INTERPOLATION_TYPE)
3734  CASE(coordinate_radial_interpolation_type)
3735  !Do nothing
3736  CASE(coordinate_radial_squared_interpolation_type)
3737  r=VALUE**(1.0_dp/3.0_dp)
3738  IF(partial_derivative_index==no_part_deriv) THEN
3739  VALUE=r
3740  ELSE
3741  VALUE=VALUE/(3.0_dp*r*r)
3742  ENDIF
3743  CASE DEFAULT
3744  local_error="The radial interpolation type of "//trim(number_to_vstring(coordinate_system% &
3745  & radial_interpolation_type,"*",err,error))//" is invalid for a cylindrical/spherical coordinate system."
3746  CALL flagerror(local_error,err,error,*999)
3747  END SELECT
3748  CASE(coordinate_prolate_spheroidal_type)
3749  SELECT CASE(coordinate_system%RADIAL_INTERPOLATION_TYPE)
3750  CASE(coordinate_radial_interpolation_type)
3751  !Do nothing
3752  CASE(coordinate_radial_squared_interpolation_type)
3753  focus=coordinate_system%FOCUS
3754  ss=VALUE/(focus*focus)
3755  sinhx=sqrt(ss)
3756  coshx=sqrt(1.0_dp+ss)
3757  IF(partial_derivative_index==no_part_deriv) THEN
3758  VALUE=log(sinhx+coshx)
3759  ELSE
3760  VALUE=VALUE/(2.0_dp*focus*focus*sinhx*coshx)
3761  ENDIF
3762  CASE(coordinate_radial_cubed_interpolation_type)
3763  focus=coordinate_system%FOCUS
3764  css=VALUE/(focus**3.0_dp)
3765  des=css*css-4.0_dp/27.0_dp
3766  IF(des>0.0_dp) THEN
3767  d=((css+sqrt(des))/2.0_dp)**(1.0_dp/3.0_dp)
3768  coshx=d+1.0_dp/(3.0_dp*d)
3769  ELSE
3770  theta=acos(css*sqrt(27.0_dp)/2.0_dp)
3771  coshx=2.0_dp/sqrt(3.0_dp)*cos(theta/3.0_dp)
3772  ENDIF
3773  sinhx=sqrt(abs(coshx*coshx-1.0_dp))
3774  IF(partial_derivative_index==no_part_deriv) THEN
3775  VALUE=log(sinhx+coshx)
3776  ELSE
3777  VALUE=VALUE/((3.0_dp*coshx*coshx-1.0_dp)*sinhx)/(focus**3.0_dp)
3778  ENDIF
3779  CASE DEFAULT
3780  local_error="The radial interpolation type of "//trim(number_to_vstring(coordinate_system% &
3781  & radial_interpolation_type,"*",err,error))//" is invalid for a prolate spheroidal coordinate system."
3782  CALL flagerror(local_error,err,error,*999)
3783  END SELECT
3784  CASE(coordinate_oblate_spheroidal_type)
3785  CALL flagerror("Not implemented.",err,error,*999)
3786  CASE DEFAULT
3787  local_error="The coordinate system type of "//trim(number_to_vstring(coordinate_system%TYPE,"*",err,error))// &
3788  & " is invalid."
3789  CALL flagerror(local_error,err,error,*999)
3790  END SELECT
3791  ELSE
3792  CALL flagerror("Coordinate system is not associated.",err,error,*999)
3793  ENDIF
3794 
3795  exits("COORDINATE_INTERPOLATION_ADJUST")
3796  RETURN
3797 999 errorsexits("COORDINATE_INTERPOLATION_ADJUST",err,error)
3798  RETURN 1
3799  END SUBROUTINE coordinate_interpolation_adjust
3800 
3801  !
3802  !================================================================================================================================
3803  !
3804 
3806  SUBROUTINE coordinate_interpolation_parameters_adjust(COORDINATE_SYSTEM,INTERPOLATION_PARAMETERS,ERR,ERROR,*)
3807 
3808  !Argument variables
3809  TYPE(coordinate_system_type), POINTER :: COORDINATE_SYSTEM
3810  TYPE(field_interpolation_parameters_type), POINTER :: INTERPOLATION_PARAMETERS
3811  INTEGER(INTG), INTENT(OUT) :: ERR
3812  TYPE(varying_string), INTENT(OUT) :: ERROR
3813  !Local variables
3814  TYPE(varying_string) :: LOCAL_ERROR
3815 
3816  enters("COORDINATE_INTERPOLATION_PARAMETERS_ADJUST",err,error,*999)
3817 
3818 !!TODO: Tidy up element parameters for non-rc coordinate systems. See bottom of XPXE and ZPZE.
3819 
3820  IF(ASSOCIATED(coordinate_system)) THEN
3821  IF(ASSOCIATED(interpolation_parameters)) THEN
3822  SELECT CASE(coordinate_system%TYPE)
3823  CASE(coordinate_rectangular_cartesian_type)
3824  !Do nothing
3825  CASE(coordinate_cylindrical_polar_type)
3826  SELECT CASE(coordinate_system%RADIAL_INTERPOLATION_TYPE)
3827  CASE(coordinate_radial_interpolation_type)
3828  !Do nothing
3829  CASE(coordinate_radial_squared_interpolation_type)
3830  CASE DEFAULT
3831  local_error="The radial interpolation type of "//trim(number_to_vstring(coordinate_system% &
3832  & radial_interpolation_type,"*",err,error))//" is invalid for a cylindrical coordinate system."
3833  CALL flagerror(local_error,err,error,*999)
3834  END SELECT
3835  CALL flagerror("Not implemented",err,error,*999)
3836  CASE(coordinate_spherical_polar_type)
3837  SELECT CASE(coordinate_system%RADIAL_INTERPOLATION_TYPE)
3838  CASE(coordinate_radial_interpolation_type)
3839  !Do nothing
3840  CASE(coordinate_radial_squared_interpolation_type)
3841  CASE DEFAULT
3842  local_error="The radial interpolation type of "//trim(number_to_vstring(coordinate_system% &
3843  & radial_interpolation_type,"*",err,error))//" is invalid for a spherical coordinate system."
3844  CALL flagerror(local_error,err,error,*999)
3845  END SELECT
3846  CALL flagerror("Not implemented.",err,error,*999)
3847  CASE(coordinate_prolate_spheroidal_type)
3848  SELECT CASE(coordinate_system%RADIAL_INTERPOLATION_TYPE)
3849  CASE(coordinate_radial_interpolation_type)
3850  !Do nothing
3851  CASE(coordinate_radial_squared_interpolation_type)
3852  CASE(coordinate_radial_cubed_interpolation_type)
3853  CASE DEFAULT
3854  local_error="The radial interpolation type of "//trim(number_to_vstring(coordinate_system% &
3855  & radial_interpolation_type,"*",err,error))//" is invalid for a prolate spheroidal coordinate system."
3856  CALL flagerror(local_error,err,error,*999)
3857  END SELECT
3858  CALL flagerror("Not implemented.",err,error,*999)
3859  CASE(coordinate_oblate_spheroidal_type)
3860  CALL flagerror("Not implemented.",err,error,*999)
3861  CASE DEFAULT
3862  local_error="The coordinate system type of "//trim(number_to_vstring(coordinate_system%TYPE,"*",err,error))// &
3863  & " is invalid."
3864  CALL flagerror(local_error,err,error,*999)
3865  END SELECT
3866  ELSE
3867  CALL flagerror("Interpolation parameters is not associated.",err,error,*999)
3868  ENDIF
3869  ELSE
3870  CALL flagerror("Coordinate system is not associated.",err,error,*999)
3871  ENDIF
3872 
3873  exits("COORDINATE_INTERPOLATION_PARAMETERS_ADJUST")
3874  RETURN
3875 999 errorsexits("COORDINATE_INTERPOLATION_PARAMETERS_ADJUST",err,error)
3876  RETURN 1
3877  END SUBROUTINE coordinate_interpolation_parameters_adjust
3878 
3879  !
3880  !================================================================================================================================
3881  !
3882 
3884  SUBROUTINE coordinates_materialsystemcalculate(geometricInterpPointMetrics,fibreInterpPoint,dNudX,dXdNu,dNudXi,dXidNu,err,error,*)
3885 
3886  !Argument variables
3887  TYPE(field_interpolated_point_metrics_type), POINTER :: geometricInterpPointMetrics
3888  TYPE(field_interpolated_point_type), POINTER :: fibreInterpPoint
3889  REAL(DP), INTENT(OUT) :: dNudX(:,:)
3890  REAL(DP), INTENT(OUT) :: dXdNu(:,:)
3891  REAL(DP), INTENT(OUT) :: dNudXi(:,:)
3892  REAL(DP), INTENT(OUT) :: dXidNu(:,:)
3893  INTEGER(INTG), INTENT(OUT) :: err
3894  TYPE(varying_string), INTENT(OUT) :: error
3895  !Local variables
3896  INTEGER(INTG) :: numberOfXDimensions,numberOfXiDimensions,numberOfNuDimensions,xiIdx
3897  REAL(DP) :: dNudXiTemp(3,3),Jnuxi
3898  TYPE(varying_string) :: localError
3899 
3900  enters("Coordinates_MaterialSystemCalculate",err,error,*999)
3901 
3902  IF(ASSOCIATED(geometricinterppointmetrics)) THEN
3903 
3904  numberofxdimensions=geometricinterppointmetrics%NUMBER_OF_X_DIMENSIONS
3905  numberofxidimensions=geometricinterppointmetrics%NUMBER_OF_XI_DIMENSIONS
3906 
3907  !Calculate dX/dNu and its inverse dNu/dX (same as transpose due to orthogonality)
3908 
3909  !The fibre interpolated point might not be used for isotropic constitutive relations
3910  IF(ASSOCIATED(fibreinterppoint)) THEN
3911  !We have a fibre field
3912  numberofnudimensions=SIZE(fibreinterppoint%values,1)
3913  SELECT CASE(numberofxdimensions)
3914  CASE(1)
3915  dxdnu(1,1)=1.0_dp
3916  CASE(2)
3917  CALL coordinates_materialsystemcalculatedxdnu2d(geometricinterppointmetrics,fibreinterppoint%values(1: &
3918  & numberofnudimensions,1),dxdnu(1:numberofxdimensions,1:numberofxdimensions),err,error,*999)
3919  CASE(3)
3920  CALL coordinates_materialsystemcalculatedxdnu3d(geometricinterppointmetrics,fibreinterppoint%values(1: &
3921  & numberofnudimensions,1),dxdnu(1:numberofxdimensions,1:numberofxdimensions),err,error,*999)
3922  CASE DEFAULT
3923  localerror="The number of dimensions in the geometric interpolated point of "// &
3924  & trim(numbertovstring(numberofxdimensions,"*",err,error))// &
3925  & " is invalid. The number of dimensions must be >= 1 and <= 3."
3926  CALL flagerror(localerror,err,error,*999)
3927  END SELECT
3928  ELSE
3929  !No fibre field
3930  numberofnudimensions=0
3931  DO xiidx=1,numberofxidimensions
3932  dnudxitemp(1:numberofxdimensions,xiidx)=geometricinterppointmetrics%DX_DXI(1:numberofxdimensions,xiidx)
3933  dxdnu(1:numberofxdimensions,xiidx)=normalise(dnudxitemp(1:numberofxdimensions,xiidx),err,error)
3934  ENDDO !xiIdx
3935  ENDIF
3936  !Calculate dNu/dX the inverse of dX/dNu (same as transpose due to orthogonality)
3937  CALL matrixtranspose(dxdnu(1:numberofxdimensions,1:numberofxdimensions),dnudx(1:numberofxdimensions,1: &
3938  & numberofxdimensions),err,error,*999)
3939  !Calculate dNu/dXi = dNu/dX * dX/dXi and its inverse dXi/dNu
3940  CALL matrixproduct(dnudx(1:numberofxdimensions,1:numberofxdimensions), &
3941  & geometricinterppointmetrics%DX_DXI(1:numberofxdimensions,1:numberofxidimensions), &
3942  & dnudxitemp(1:numberofxdimensions,1:numberofxidimensions),err,error,*999)
3943  !Setup dNudXi
3944  CALL identitymatrix(dnudxi,err,error,*999)
3945  dnudxi(1:numberofxdimensions,1:numberofxidimensions)=dnudxitemp(1:numberofxdimensions,1:numberofxidimensions)
3946 
3947  IF(numberofxdimensions==numberofxidimensions) THEN
3948  CALL invert(dnudxi(1:numberofxdimensions,1:numberofxidimensions),dxidnu(1:numberofxidimensions,1:numberofxdimensions), &
3949  & jnuxi,err,error,*999)
3950  ELSE
3951  CALL flagerror("Not implemented",err,error,*999)
3952  ENDIF
3953 
3954  IF(diagnostics1) THEN
3955  CALL writestring(diagnostic_output_type,"",err,error,*999)
3956  CALL writestring(diagnostic_output_type,"Calculated material coordinate system:",err,error,*999)
3957  CALL writestringvalue(diagnostic_output_type," Number of X dimensions = ",numberofxdimensions,err,error,*999)
3958  CALL writestringvalue(diagnostic_output_type," Number of Xi dimensions = ",numberofxidimensions,err,error,*999)
3959  CALL writestringvalue(diagnostic_output_type," Number of Nu dimensions = ",numberofnudimensions,err,error,*999)
3960  CALL writestring(diagnostic_output_type," Derivative of X wrt Nu:",err,error,*999)
3961  CALL writestringmatrix(diagnostic_output_type,1,1,numberofxdimensions,1,1,numberofxdimensions, &
3962  & numberofxdimensions,numberofxdimensions,dxdnu,write_string_matrix_name_and_indices, &
3963  & '(" dX_dNu','(",I1,",:)',' :",3(X,E13.6))','(18X,3(X,E13.6))',err,error,*999)
3964  CALL writestring(diagnostic_output_type," Derivative of Nu wrt X:",err,error,*999)
3965  CALL writestringmatrix(diagnostic_output_type,1,1,numberofxdimensions,1,1,numberofxdimensions, &
3966  & numberofxdimensions,numberofxdimensions,dnudx,write_string_matrix_name_and_indices, &
3967  & '(" dNu_dX','(",I1,",:)',' :",3(X,E13.6))','(18X,3(X,E13.6))',err,error,*999)
3968  CALL writestringvalue(diagnostic_output_type," Determinant dNu_dX, JNuX = ", &
3969  & determinant(dnudx(1:numberofxdimensions,1:numberofxdimensions),err,error),err,error,*999)
3970  CALL writestring(diagnostic_output_type," Derivative of Nu wrt Xi:",err,error,*999)
3971  CALL writestringmatrix(diagnostic_output_type,1,1,numberofxdimensions,1,1,numberofxidimensions, &
3972  & numberofxidimensions,numberofxidimensions,dnudxi,write_string_matrix_name_and_indices, &
3973  & '(" dNu_dXi','(",I1,",:)',' :",3(X,E13.6))','(18X,3(X,E13.6))',err,error,*999)
3974  CALL writestring(diagnostic_output_type," Derivative of Xi wrt Nu:",err,error,*999)
3975  CALL writestringmatrix(diagnostic_output_type,1,1,numberofxidimensions,1,1,numberofxdimensions, &
3976  & numberofxdimensions,numberofxdimensions,dxidnu,write_string_matrix_name_and_indices, &
3977  & '(" dXi_dNu','(",I1,",:)',' :",3(X,E13.6))','(18X,3(X,E13.6))',err,error,*999)
3978  CALL writestringvalue(diagnostic_output_type," Determinant dNu_dXi, JNuXi = ",jnuxi,err,error,*999)
3979  ENDIF
3980 
3981  ELSE
3982  CALL flagerror("Geometric interpolated point metrics is not associated.",err,error,*999)
3983  ENDIF
3984 
3985  exits("Coordinates_MaterialSystemCalculate")
3986  RETURN
3987 999 errorsexits("Coordinates_MaterialSystemCalculate",err,error)
3988  RETURN 1
3989 
3990  END SUBROUTINE coordinates_materialsystemcalculate
3991 
3992  !
3993  !================================================================================================================================
3994  !
3995 
3997  SUBROUTINE coordinates_materialsystemcalculatedxdnu2d(geometricInterpPointMetrics,angle,dXdNu,err,error,*)
3998 
3999  !Argument variables
4000  TYPE(field_interpolated_point_metrics_type), POINTER :: geometricInterpPointMetrics
4001  REAL(DP), INTENT(IN) :: angle(:)
4002  REAL(DP), INTENT(OUT) :: dXdNu(:,:)
4003  INTEGER(INTG), INTENT(OUT) :: err
4004  TYPE(varying_string), INTENT(OUT) :: error
4005  !Local Variables
4006  REAL(DP) :: det,dXdNuR(2,2),R(2,2),f(2),g(2)
4007 
4008  enters("Coordinates_MaterialSystemCalculatedXdNu2D",err,error,*999)
4009 
4010  IF(ASSOCIATED(geometricinterppointmetrics)) THEN
4011 
4012  !First calculate reference material CS
4013 
4014  !Reference material direction 1.
4015  f(1:2) = [ geometricinterppointmetrics%DX_DXI(1,1),geometricinterppointmetrics%DX_DXI(2,1) ]
4016 
4017  !Compute (normalised) vector orthogonal to material direction 1 to form material direction 2
4018  g(1:2) = [ -1.0_dp*f(2),f(1) ]
4019 
4020  dxdnur(1:2,1) = normalise(f(1:2),err,error)
4021  IF(err/=0) GOTO 999
4022  dxdnur(1:2,2) = normalise(g(1:2),err,error)
4023  IF(err/=0) GOTO 999
4024 
4025  !Rotate by multiply with rotation matrix
4026  r(:,1) = [ cos(angle(1)),-1.0_dp*sin(angle(1)) ]
4027  r(:,2) = [ sin(angle(1)),cos(angle(1)) ]
4028 
4029  CALL matrixproduct(r,dxdnur,dxdnu,err,error,*999)
4030 
4031  dxdnu(1:2,1) = normalise(dxdnu(1:2,1),err,error)
4032  IF(err/=0) GOTO 999
4033  dxdnu(1:2,2) = normalise(dxdnu(1:2,2),err,error)
4034  IF(err/=0) GOTO 999
4035 
4036  IF(diagnostics1) THEN
4037  CALL writestring(diagnostic_output_type,"",err,error,*999)
4038  CALL writestring(diagnostic_output_type,"2D material system calculation:",err,error,*999)
4039  CALL writestring(diagnostic_output_type," Reference material directions:",err,error,*999)
4040  CALL writestringvector(diagnostic_output_type,1,1,2,2,2,f,'(" f :",2(X,E13.6))','(20X,2(X,E13.6))', &
4041  & err,error,*999)
4042  CALL writestringvector(diagnostic_output_type,1,1,2,2,2,g,'(" g :",2(X,E13.6))','(20X,2(X,E13.6))', &
4043  & err,error,*999)
4044  CALL writestring(diagnostic_output_type," Derivative of X wrt Nu (reference):",err,error,*999)
4045  CALL writestringmatrix(diagnostic_output_type,1,1,2,1,1,2,2,2,dxdnur,write_string_matrix_name_and_indices, &
4046  & '(" dX_dNuR','(",I1,",:)',' :",2(X,E13.6))','(20X,2(X,E13.6))',err,error,*999)
4047  CALL writestring(diagnostic_output_type," Fibre calculation:",err,error,*999)
4048  CALL writestringvalue(diagnostic_output_type," Fibre angle = ",angle(1),err,error,*999)
4049  IF(diagnostics2) THEN
4050  CALL writestring(diagnostic_output_type," Rotation matrix:",err,error,*999)
4051  CALL writestringmatrix(diagnostic_output_type,1,1,2,1,1,2,2,2,r,write_string_matrix_name_and_indices, &
4052  & '(" R','(",I1,",:)',' :",2(X,E13.6))','(20X,2(X,E13.6))',err,error,*999)
4053  ENDIF
4054  CALL writestring(diagnostic_output_type," Derivative of X wrt Nu (material):",err,error,*999)
4055  CALL writestringmatrix(diagnostic_output_type,1,1,2,1,1,2,2,2,dxdnu,write_string_matrix_name_and_indices, &
4056  & '(" dX_dNu','(",I1,",:)',' :",2(X,E13.6))','(20X,2(X,E13.6))',err,error,*999)
4057  det=determinant(dxdnu,err,error)
4058  CALL writestringvalue(diagnostic_output_type," Determinant dX_dNu = ",det,err,error,*999)
4059  ENDIF
4060 
4061  ELSE
4062  CALL flagerror("Geometry interpolated point metrics is not associated.",err,error,*999)
4063  ENDIF
4064 
4065  exits("Coordinates_MaterialSystemCalculatedXdNu2D")
4066  RETURN
4067 999 errorsexits("Coordinates_MaterialSystemCalculatedXdNu2D",err,error)
4068  RETURN 1
4069 
4070  END SUBROUTINE coordinates_materialsystemcalculatedxdnu2d
4071 
4072  !
4073  !================================================================================================================================
4074  !
4075 
4077  SUBROUTINE coordinates_materialsystemcalculatedxdnu3d(geometricInterpPointMetrics,angle,dXdNu,err,error,*)
4078 
4079  !Argument variables
4080  TYPE(field_interpolated_point_metrics_type), POINTER :: geometricInterpPointMetrics
4081  REAL(DP), INTENT(IN) :: angle(:)
4082  REAL(DP), INTENT(OUT) :: dXdNu(:,:)
4083  INTEGER(INTG), INTENT(OUT) :: err
4084  TYPE(varying_string), INTENT(OUT) :: error
4085  !Local Variables
4086  REAL(DP) :: angles(3),det,dXdNu2(3,3),dXdNu3(3,3),dXdNuR(3,3),f(3),g(3),h(3), &
4087  & Raf(3,3),Rbf(3,3),Rai(3,3),Rbi(3,3),Ras(3,3),Rbs(3,3)
4088 
4089  enters("Coordinates_MaterialSystemCalculatedXdNu3D",err,error,*999)
4090 
4091  IF(ASSOCIATED(geometricinterppointmetrics)) THEN
4092 
4093  !Calculate fibre, imbrication and sheet angles (allowing for missing angles)
4094  angles=0.0_dp
4095  angles(1:SIZE(angle,1))=angle(1:SIZE(angle,1))
4096 
4097  !First calculate reference material CS
4098  f(1:3)=geometricinterppointmetrics%DX_DXI(1:3,1) !reference material direction 1.
4099  CALL crossproduct(geometricinterppointmetrics%DX_DXI(1:3,1),geometricinterppointmetrics%DX_DXI(1:3,2),h,err,error,*999) !reference material direction 3.
4100  CALL crossproduct(h,f,g,err,error,*999) !reference material direction 2.
4101 
4102  dxdnur(1:3,1)=normalise(f,err,error)
4103  IF(err/=0) GOTO 999
4104  dxdnur(1:3,2)=normalise(g,err,error)
4105  IF(err/=0) GOTO 999
4106  dxdnur(1:3,3)=normalise(h,err,error)
4107  IF(err/=0) GOTO 999
4108 
4109  IF(diagnostics1) THEN
4110  ENDIF
4111 
4112  !FIBRE ANGLE(alpha) - angles(1)
4113  !In order to rotate reference material CS by alpha(fibre angle) in anti-clockwise
4114  !direction about its axis 3, following steps are performed.
4115  !(a) first align reference material direction 3 with Z(spatial) axis by rotating the ref material CS.
4116  !(b) then rotate the aligned material CS by alpha about Z axis in anti-clockwise direction
4117  !(c) apply the inverse of step(a) to the CS in (b)
4118  !It can be shown that steps (a),(b) and (c) are equivalent to post-multiplying
4119  !rotation in (a) by rotation in (b). i.e. Ra*Rb
4120 
4121  !The normalised reference material CS contains the transformation(rotation) between
4122  !the spatial CS -> reference material CS. i.e. Raf
4123  raf=dxdnur
4124 
4125  !Initialise rotation matrix Rbf
4126  CALL identitymatrix(rbf,err,error,*999)
4127  !Populate rotation matrix Rbf about axis 3 (Z)
4128  rbf(1,1)=cos(angles(1))
4129  rbf(1,2)=-1.0_dp*sin(angles(1))
4130  rbf(2,1)=sin(angles(1))
4131  rbf(2,2)=cos(angles(1))
4132 
4133  CALL matrixproduct(raf,rbf,dxdnu3,err,error,*999)
4134 
4135  !IMBRICATION ANGLE (beta) - angles(2)
4136  !In order to rotate alpha-rotated material CS by beta(imbrication angle) in anti-clockwise
4137  !direction about its new axis 2, following steps are performed.
4138  !(a) first align new material direction 2 with Y(spatial) axis by rotating the new material CS.
4139  !(b) then rotate the aligned CS by beta about Y axis in anti-clockwise direction
4140  !(c) apply the inverse of step(a) to the CS in (b)
4141  !As mentioned above, (a),(b) and (c) are equivalent to post-multiplying
4142  !rotation in (a) by rotation in (b). i.e. Rai*Rbi
4143 
4144  !dXdNu3 contains the transformation(rotation) between
4145  !the spatial CS -> alpha-rotated reference material CS. i.e. Rai
4146  rai=dxdnu3
4147  !Initialise rotation matrix Rbi
4148  CALL identitymatrix(rbi,err,error,*999)
4149  !Populate rotation matrix Rbi about axis 2 (Y). Note the sign change
4150  rbi(1,1)=cos(angles(2))
4151  rbi(1,3)=sin(angles(2))
4152  rbi(3,1)=-1.0_dp*sin(angles(2))
4153  rbi(3,3)=cos(angles(2))
4154 
4155  CALL matrixproduct(rai,rbi,dxdnu2,err,error,*999)
4156 
4157  !SHEET ANGLE (gamma) - angles(3)
4158  !In order to rotate alpha-beta-rotated material CS by gama(sheet angle) in anti-clockwise
4159  !direction about its new axis 1, following steps are performed.
4160  !(a) first align new material direction 1 with X(spatial) axis by rotating the new material CS.
4161  !(b) then rotate the aligned CS by gama about X axis in anti-clockwise direction
4162  !(c) apply the inverse of step(a) to the CS in (b)
4163  !Again steps (a),(b) and (c) are equivalent to post-multiplying
4164  !rotation in (a) by rotation in (b). i.e. Ras*Rbs
4165 
4166  !dXdNu2 contains the transformation(rotation) between
4167  !the spatial CS -> alpha-beta-rotated reference material CS. i.e. Ras
4168  ras=dxdnu2
4169  !Initialise rotation matrix Rbs
4170  CALL identitymatrix(rbs,err,error,*999)
4171  !Populate rotation matrix Rbs about axis 1 (X).
4172  rbs(2,2)=cos(angles(3))
4173  rbs(2,3)=-1.0_dp*sin(angles(3))
4174  rbs(3,2)=sin(angles(3))
4175  rbs(3,3)=cos(angles(3))
4176 
4177  CALL matrixproduct(ras,rbs,dxdnu,err,error,*999)
4178 
4179  IF(diagnostics1) THEN
4180  CALL writestring(diagnostic_output_type,"",err,error,*999)
4181  CALL writestring(diagnostic_output_type,"3D material system calculation:",err,error,*999)
4182  CALL writestring(diagnostic_output_type," Reference material directions:",err,error,*999)
4183  CALL writestringvector(diagnostic_output_type,1,1,3,3,3,f,'(" f :",3(X,E13.6))','(20X,3(X,E13.6))', &
4184  & err,error,*999)
4185  CALL writestringvector(diagnostic_output_type,1,1,3,3,3,g,'(" g :",3(X,E13.6))','(20X,3(X,E13.6))', &
4186  & err,error,*999)
4187  CALL writestringvector(diagnostic_output_type,1,1,3,3,3,h,'(" h :",3(X,E13.6))','(20X,3(X,E13.6))', &
4188  & err,error,*999)
4189  CALL writestring(diagnostic_output_type," Derivative of X wrt Nu (reference):",err,error,*999)
4190  CALL writestringmatrix(diagnostic_output_type,1,1,3,1,1,3,3,3,dxdnur,write_string_matrix_name_and_indices, &
4191  & '(" dX_dNuR','(",I1,",:)',' :",3(X,E13.6))','(20X,3(X,E13.6))',err,error,*999)
4192  CALL writestring(diagnostic_output_type," Fibre calculation:",err,error,*999)
4193  CALL writestringvalue(diagnostic_output_type," Fibre angle = ",angles(1),err,error,*999)
4194  IF(diagnostics2) THEN
4195  CALL writestring(diagnostic_output_type," Rotation matrix A:",err,error,*999)
4196  CALL writestringmatrix(diagnostic_output_type,1,1,3,1,1,3,3,3,raf,write_string_matrix_name_and_indices, &
4197  & '(" Ra','(",I1,",:)',' :",3(X,E13.6))','(20X,3(X,E13.6))',err,error,*999)
4198  CALL writestring(diagnostic_output_type," Rotation matrix B:",err,error,*999)
4199  CALL writestringmatrix(diagnostic_output_type,1,1,3,1,1,3,3,3,rbf,write_string_matrix_name_and_indices, &
4200  & '(" Rb','(",I1,",:)',' :",3(X,E13.6))','(20X,3(X,E13.6))',err,error,*999)
4201  ENDIF
4202  CALL writestring(diagnostic_output_type," Derivative of X wrt Nu (after alpha rotation):",err,error,*999)
4203  CALL writestringmatrix(diagnostic_output_type,1,1,3,1,1,3,3,3,dxdnu3,write_string_matrix_name_and_indices, &
4204  & '(" dX_dNu3','(",I1,",:)',' :",3(X,E13.6))','(20X,3(X,E13.6))',err,error,*999)
4205  CALL writestring(diagnostic_output_type," Imbrication calculation:",err,error,*999)
4206  CALL writestringvalue(diagnostic_output_type," Imbrication angle = ",angles(2),err,error,*999)
4207  IF(diagnostics2) THEN
4208  CALL writestring(diagnostic_output_type," Rotation matrix A:",err,error,*999)
4209  CALL writestringmatrix(diagnostic_output_type,1,1,3,1,1,3,3,3,rai,write_string_matrix_name_and_indices, &
4210  & '(" Ra','(",I1,",:)',' :",3(X,E13.6))','(20X,3(X,E13.6))',err,error,*999)
4211  CALL writestring(diagnostic_output_type," Rotation matrix B:",err,error,*999)
4212  CALL writestringmatrix(diagnostic_output_type,1,1,3,1,1,3,3,3,rbi,write_string_matrix_name_and_indices, &
4213  & '(" Rb','(",I1,",:)',' :",3(X,E13.6))','(20X,3(X,E13.6))',err,error,*999)
4214  ENDIF
4215  CALL writestring(diagnostic_output_type," Derivative of X wrt Nu (after alpha-beta rotation):",err,error,*999)
4216  CALL writestringmatrix(diagnostic_output_type,1,1,3,1,1,3,3,3,dxdnu2,write_string_matrix_name_and_indices, &
4217  & '(" dX_dNu2','(",I1,",:)',' :",3(X,E13.6))','(20X,3(X,E13.6))',err,error,*999)
4218  CALL writestring(diagnostic_output_type," Sheet calculation:",err,error,*999)
4219  CALL writestringvalue(diagnostic_output_type," Sheet angle = ",angles(3),err,error,*999)
4220  IF(diagnostics2) THEN
4221  CALL writestring(diagnostic_output_type," Rotation matrix A:",err,error,*999)
4222  CALL writestringmatrix(diagnostic_output_type,1,1,3,1,1,3,3,3,ras,write_string_matrix_name_and_indices, &
4223  & '(" Ra','(",I1,",:)',' :",3(X,E13.6))','(20X,3(X,E13.6))',err,error,*999)
4224  CALL writestring(diagnostic_output_type," Rotation matrix B:",err,error,*999)
4225  CALL writestringmatrix(diagnostic_output_type,1,1,3,1,1,3,3,3,rbs,write_string_matrix_name_and_indices, &
4226  & '(" Rb','(",I1,",:)',' :",3(X,E13.6))','(20X,3(X,E13.6))',err,error,*999)
4227  ENDIF
4228  CALL writestring(diagnostic_output_type," Derivative of X wrt Nu (material):",err,error,*999)
4229  CALL writestringmatrix(diagnostic_output_type,1,1,3,1,1,3,3,3,dxdnu,write_string_matrix_name_and_indices, &
4230  & '(" dX_dNu','(",I1,",:)',' :",3(X,E13.6))','(20X,3(X,E13.6))',err,error,*999)
4231  det=determinant(dxdnu,err,error)
4232  CALL writestringvalue(diagnostic_output_type," Determinant dX_dNu = ",det,err,error,*999)
4233  ENDIF
4234 
4235  ELSE
4236  CALL flagerror("Geometry interpolated point metrics is not associated.",err,error,*999)
4237  ENDIF
4238 
4239  exits("Coordinates_MaterialSystemCalculatedXdNu3D")
4240  RETURN
4241 999 errorsexits("Coordinates_MaterialSystemCalculatedXdNu3D",err,error)
4242  RETURN 1
4243 
4244  END SUBROUTINE coordinates_materialsystemcalculatedxdnu3d
4245 
4246  !
4247  !================================================================================================================================
4248  !
4249 
4252  SUBROUTINE coordinate_system_user_number_find(USER_NUMBER,COORDINATE_SYSTEM,ERR,ERROR,*)
4253 
4254  !Argument variables
4255  INTEGER(INTG), INTENT(IN) :: USER_NUMBER
4256  TYPE(coordinate_system_type), POINTER :: COORDINATE_SYSTEM
4257  INTEGER(INTG), INTENT(OUT) :: ERR
4258  TYPE(varying_string), INTENT(OUT) :: ERROR
4259  !Local Variables
4260  INTEGER(INTG) :: coord_system_idx
4261 
4262  enters("COORDINATE_SYSTEM_USER_NUMBER_FIND",err,error,*999)
4263 
4264  IF(ASSOCIATED(coordinate_system)) THEN
4265  CALL flagerror("Coordinate_system is already associated.",err,error,*999)
4266  ELSE
4267  NULLIFY(coordinate_system)
4268  coord_system_idx=1
4269  DO WHILE(coord_system_idx<=coordinate_systems%NUMBER_OF_COORDINATE_SYSTEMS.AND..NOT.ASSOCIATED(coordinate_system))
4270  IF(coordinate_systems%COORDINATE_SYSTEMS(coord_system_idx)%PTR%USER_NUMBER==user_number) THEN
4271  coordinate_system=>coordinate_systems%COORDINATE_SYSTEMS(coord_system_idx)%PTR
4272  ELSE
4273  coord_system_idx=coord_system_idx+1
4274  ENDIF
4275  ENDDO
4276  ENDIF
4277 
4278  exits("COORDINATE_SYSTEM_USER_NUMBER_FIND")
4279  RETURN
4280 999 errorsexits("COORDINATE_SYSTEM_USER_NUMBER_FIND",err,error)
4281  RETURN 1
4282  END SUBROUTINE coordinate_system_user_number_find
4283 
4284  !
4285  !================================================================================================================================
4286  !
4287 
4289  SUBROUTINE coordinate_systems_finalise(ERR,ERROR,*)
4290 
4291  !Argument variables
4292  INTEGER(INTG), INTENT(OUT) :: ERR
4293  TYPE(varying_string), INTENT(OUT) :: ERROR
4294  !Local Variables
4295  INTEGER(INTG) :: coord_system_idx
4296 
4297  enters("COORDINATE_SYSTEMS_FINALISE",err,error,*999)
4298 
4299  DO coord_system_idx=1,coordinate_systems%NUMBER_OF_COORDINATE_SYSTEMS
4300  CALL coordinate_system_finalise(coordinate_systems%COORDINATE_SYSTEMS(coord_system_idx)%PTR,err,error,*999)
4301  ENDDO !coord_system_idx
4302  DEALLOCATE(coordinate_systems%COORDINATE_SYSTEMS)
4303  coordinate_systems%NUMBER_OF_COORDINATE_SYSTEMS=0
4304 
4305  exits("COORDINATE_SYSTEMS_FINALISE")
4306  RETURN
4307 999 errorsexits("COORDINATE_SYSTEMS_FINALISE",err,error)
4308  RETURN 1
4309  END SUBROUTINE coordinate_systems_finalise
4310 
4311  !
4312  !================================================================================================================================
4313  !
4314 
4316  SUBROUTINE coordinate_systems_initialise(ERR,ERROR,*)
4317 
4318  !Argument variables
4319  INTEGER(INTG), INTENT(OUT) :: ERR
4320  TYPE(varying_string), INTENT(OUT) :: ERROR
4321  !Local Variables
4322 
4323  enters("COORDINATE_SYSTEMS_INITIALISE",err,error,*999)
4324 
4325  !Allocate the coordinate systems
4326  ALLOCATE(coordinate_systems%COORDINATE_SYSTEMS(1),stat=err)
4327  IF(err/=0) CALL flagerror("Could not allocate coordinate systems.",err,error,*999)
4328  !Create the default RC World cooordinate system
4329  ALLOCATE(coordinate_systems%COORDINATE_SYSTEMS(1)%PTR,stat=err)
4330  IF(err/=0) CALL flagerror("Could not allocate world coordinate system.",err,error,*999)
4331  coordinate_systems%COORDINATE_SYSTEMS(1)%PTR%USER_NUMBER=0
4332  coordinate_systems%COORDINATE_SYSTEMS(1)%PTR%TYPE=coordinate_rectangular_cartesian_type
4333  coordinate_systems%COORDINATE_SYSTEMS(1)%PTR%NUMBER_OF_DIMENSIONS=3
4334  coordinate_systems%COORDINATE_SYSTEMS(1)%PTR%FOCUS=1.0_dp
4335  coordinate_systems%COORDINATE_SYSTEMS(1)%PTR%ORIGIN=(/0.0_dp,0.0_dp,0.0_dp/)
4336  coordinate_systems%COORDINATE_SYSTEMS(1)%PTR%ORIENTATION=reshape(&
4337  & (/1.0_dp,0.0_dp,0.0_dp, &
4338  & 0.0_dp,1.0_dp,0.0_dp, &
4339  & 0.0_dp,0.0_dp,1.0_dp/), &
4340  & (/3,3/))
4341  coordinate_systems%COORDINATE_SYSTEMS(1)%PTR%COORDINATE_SYSTEM_FINISHED=.true.
4342  coordinate_systems%NUMBER_OF_COORDINATE_SYSTEMS=1
4343 
4344  exits("COORDINATE_SYSTEMS_INITIALISE")
4345  RETURN
4346 999 errorsexits("COORDINATE_SYSTEMS_INITIALISE",err,error)
4347  RETURN 1
4348  END SUBROUTINE coordinate_systems_initialise
4349 
4350  !
4351  !================================================================================================================================
4352  !
4353 
4354 END MODULE coordinate_routines
4355 
base_routines::enters
subroutine, public enters(NAME, ERR, ERROR,)
Records the entry into the named procedure and initialises the error code.
Definition: base_routines.f90:363
input_output::write_string_value
Write a string followed by a value to a given output stream.
Definition: input_output.f90:86
coordinate_routines::coordinate_convert_to_rc_dp
real(dp) function, dimension(size(x, 1)) coordinate_convert_to_rc_dp(COORDINATE_SYSTEM, X, ERR, ERROR)
COORDINATE_CONVERT_TO_RC_DP performs a coordinate transformation from a coordinate system identified ...
Definition: coordinate_routines.f90:445
coordinate_routines::coordinate_no_radial_interpolation_type
integer(intg), parameter, public coordinate_no_radial_interpolation_type
No radial interpolation.
Definition: coordinate_routines.f90:81
coordinate_routines
This module contains all coordinate transformation and support routines.
Definition: coordinate_routines.f90:45
coordinate_routines::coordinate_prolate_spheroidal_type
integer(intg), parameter, public coordinate_prolate_spheroidal_type
Prolate spheroidal coordinate system type.
Definition: coordinate_routines.f90:73
coordinate_routines::dzx
Definition: coordinate_routines.f90:151
maths::invert
Returns the inverse of a matrix.
Definition: maths.f90:131
coordinate_routines::coordinate_derivative_convert_to_rc_sp
subroutine coordinate_derivative_convert_to_rc_sp(COORDINATE_SYSTEM, PART_DERIV_TYPE, X, Z, ERR, ERROR,)
Definition: coordinate_routines.f90:3037
constants::zero_tolerance_sp
real(sp), parameter zero_tolerance_sp
The zero tolerance for single precision zero tests i.e., if(abs(x)>zero_tolerance) then...
Definition: constants.f90:75
strings::number_to_vstring
Converts a number to its equivalent varying string representation.
Definition: strings.f90:161
coordinate_routines::coordinate_system_origin_get
subroutine, public coordinate_system_origin_get(COORDINATE_SYSTEM, ORIGIN, ERR, ERROR,)
Returns the origin of a coordinate system.
Definition: coordinate_routines.f90:1584
coordinate_routines::d2zx_dp
real(dp) function, dimension(size(x, 1)) d2zx_dp(COORDINATE_SYSTEM, I, J, X, ERR, ERROR)
Definition: coordinate_routines.f90:2038
constants::pi
real(dp), parameter pi
The double precision value of pi.
Definition: constants.f90:57
constants::no_part_deriv
integer(intg), parameter no_part_deriv
No partial derivative i.e., u.
Definition: constants.f90:177
coordinate_routines::coordinate_jacobian_area_type
integer(intg), parameter, public coordinate_jacobian_area_type
Area type Jacobian.
Definition: coordinate_routines.f90:93
coordinate_routines::coordinate_systems_initialise
subroutine, public coordinate_systems_initialise(ERR, ERROR,)
Initialises the coordinate systems and creates the world coordinate system.
Definition: coordinate_routines.f90:4317
input_output::writestringmatrix
Write a string followed by a matrix to a specified output stream.
Definition: input_output.f90:335
constants::partial_derivative_second_derivative_map
integer(intg), dimension(4) partial_derivative_second_derivative_map
PARTIAL_DERIVATIVE_SECOND_DERIVATIVE_MAP(nic) gives the partial derivative index for the second deriv...
Definition: constants.f90:256
coordinate_routines::coordinate_interpolation_parameters_adjust
subroutine, public coordinate_interpolation_parameters_adjust(COORDINATE_SYSTEM, INTERPOLATION_PARAMETERS, ERR, ERROR,)
Adjusts the interpolation parameters for non-rectangular cartesian coordinate systems.
Definition: coordinate_routines.f90:3807
strings
This module contains all string manipulation and transformation routines.
Definition: strings.f90:45
coordinate_routines::coordinate_jacobian_line_type
integer(intg), parameter, public coordinate_jacobian_line_type
Line type Jacobian.
Definition: coordinate_routines.f90:92
constants::first_part_deriv
integer(intg), parameter first_part_deriv
First partial derivative i.e., du/ds.
Definition: constants.f90:178
types::coordinate_system_ptr_type
A buffer type to allow for an array of pointers to a COORDINATE_SYSTEM_TYPE.
Definition: types.f90:267
coordinate_routines::coordinate_system_normal_calculate
subroutine coordinate_system_normal_calculate(COORDINATE_SYSTEM, REVERSE, X, N, ERR, ERROR,)
Calculates the normal vector, N, at the point X. IF REVERSE is true the reversed normal is returned...
Definition: coordinate_routines.f90:1044
constants::part_deriv_s2
integer(intg), parameter part_deriv_s2
First partial derivative in the s2 direction i.e., du/ds2.
Definition: constants.f90:183
coordinate_routines::coordinate_system_finalise
subroutine coordinate_system_finalise(COORDINATE_SYSTEM, ERR, ERROR,)
Finalises a coordinate system and deallocates all memory.
Definition: coordinate_routines.f90:1220
coordinate_routines::coordinate_derivative_norm
subroutine, public coordinate_derivative_norm(COORDINATE_SYSTEM, PART_DERIV_INDEX, INTERPOLATED_POINT, DERIV_NORM, ERR, ERROR,)
Calculates the norm of a derivative in a coordinate system identified by COORDINATE_SYSTEM at the giv...
Definition: coordinate_routines.f90:3607
coordinate_routines::coordinate_jacobian_no_type
integer(intg), parameter, public coordinate_jacobian_no_type
No Jacobian.
Definition: coordinate_routines.f90:91
coordinate_routines::coordinate_convert_from_rc_sp
real(sp) function, dimension(size(z, 1)) coordinate_convert_from_rc_sp(COORDINATE_SYSTEM, Z, ERR, ERROR)
COORDINATE_CONVERT_FROM_RC_SP performs a coordinate transformation from a rectangular cartesian coord...
Definition: coordinate_routines.f90:325
maths::normalise
Normalises a vector.
Definition: maths.f90:221
maths
This module contains all mathematics support routines.
Definition: maths.f90:45
input_output::write_string_matrix
Write a string followed by a matrix to a specified output stream.
Definition: input_output.f90:326
base_routines::diagnostics2
logical, save, public diagnostics2
.TRUE. if level 2 diagnostic output is active in the current routine
Definition: base_routines.f90:154
coordinate_routines::coordinate_system_dimension_set
subroutine, public coordinate_system_dimension_set(COORDINATE_SYSTEM, DIMENSION, ERR, ERROR,)
Sets/changes the dimension of the coordinate system.
Definition: coordinate_routines.f90:1356
coordinate_routines::coordinate_system_user_number_find
subroutine, public coordinate_system_user_number_find(USER_NUMBER, COORDINATE_SYSTEM, ERR, ERROR,)
Returns a pointer to the coordinate system identified by USER_NUMBER. If a coordinate system with tha...
Definition: coordinate_routines.f90:4253
coordinate_routines::coordinate_systems
type(coordinate_systems_type) coordinate_systems
Definition: coordinate_routines.f90:113
coordinate_routines::coordinate_system_create_start
subroutine, public coordinate_system_create_start(USER_NUMBER, COORDINATE_SYSTEM, ERR, ERROR,)
Starts the creation of and initialises a new coordinate system.
Definition: coordinate_routines.f90:1737
iso_varying_string
This module provides an iso_varying_string module, conformant to the API specified in ISO/IEC 1539-2:...
Definition: iso_varying_string.f90:62
types::coordinate_system_type
Contains information on a coordinate system.
Definition: types.f90:255
iso_varying_string::varying_string
Definition: iso_varying_string.f90:75
constants
This module contains all program wide constants.
Definition: constants.f90:45
maths::norm_cross_product
Calculates the normalised vector cross product of two vectors.
Definition: maths.f90:227
constants::part_deriv_s1
integer(intg), parameter part_deriv_s1
First partial derivative in the s1 direction i.e., du/ds1.
Definition: constants.f90:181
base_routines::flag_warning
Flags a warning to the user.
Definition: base_routines.f90:232
coordinate_routines::coordinate_system_origin_set
subroutine, public coordinate_system_origin_set(COORDINATE_SYSTEM, ORIGIN, ERR, ERROR,)
Sets/changes the origin of a coordinate system.
Definition: coordinate_routines.f90:1621
types::field_interpolated_point_metrics_type
Contains the interpolated point coordinate metrics. Old CMISS name GL,GU,RG.
Definition: types.f90:1112
coordinate_routines::coordinate_system_focus_get
subroutine, public coordinate_system_focus_get(COORDINATE_SYSTEM, FOCUS, ERR, ERROR,)
Returns the coordinate system focus.
Definition: coordinate_routines.f90:1246
coordinate_routines::coordinate_system_dimension_get
subroutine, public coordinate_system_dimension_get(COORDINATE_SYSTEM, NUMBER_OF_DIMENSIONS, ERR, ERROR,)
Gets the coordinate system dimension.
Definition: coordinate_routines.f90:1187
coordinate_routines::dxz_dp
real(dp) function, dimension(size(x, 1)) dxz_dp(COORDINATE_SYSTEM, I, X, ERR, ERROR)
Calculates DX(:)/DZ(I) at X, where Z(I) are rectangular cartesian and X(:) are curvilinear coordinate...
Definition: coordinate_routines.f90:1906
coordinate_routines::coordinate_derivative_convert_to_rc_dp
subroutine coordinate_derivative_convert_to_rc_dp(COORDINATE_SYSTEM, PART_DERIV_TYPE, X, Z, ERR, ERROR,)
Definition: coordinate_routines.f90:2469
coordinate_routines::coordinates_materialsystemcalculate
subroutine, public coordinates_materialsystemcalculate(geometricInterpPointMetrics, fibreInterpPoint, dNudX, dXdNu, dNudXi, dXidNu, err, error,)
Calculates the tensor to get from material coordinate system, nu, to local coordinate system...
Definition: coordinate_routines.f90:3885
coordinate_routines::coordinate_radial_cubed_interpolation_type
integer(intg), parameter, public coordinate_radial_cubed_interpolation_type
r^3 radial interpolation
Definition: coordinate_routines.f90:84
coordinate_routines::coordinates_radialinterpolationtypeget
subroutine, public coordinates_radialinterpolationtypeget(coordinateSystem, radialInterpolationType, err, error,)
Gets the coordinate system radial interpolation type.
Definition: coordinate_routines.f90:1285
coordinate_routines::coordinate_interpolation_adjust
subroutine, public coordinate_interpolation_adjust(COORDINATE_SYSTEM, PARTIAL_DERIVATIVE_INDEX, VALUE, ERR, ERROR,)
Adjusts the interpolation for non-rectangular cartesian coordinate systems.
Definition: coordinate_routines.f90:3699
base_routines::exits
subroutine, public exits(NAME)
Records the exit out of the named procedure.
Definition: base_routines.f90:504
coordinate_routines::coordinate_system_orientation_get
subroutine, public coordinate_system_orientation_get(COORDINATE_SYSTEM, ORIENTATION, ERR, ERROR,)
Returns the orientation of a coordinate system.
Definition: coordinate_routines.f90:1657
types
This module contains all type definitions in order to avoid cyclic module references.
Definition: types.f90:70
constants::part_deriv_s3
integer(intg), parameter part_deriv_s3
First partial derivative in the s3 direction i.e., du/ds3.
Definition: constants.f90:186
input_output::write_string
Write a string to a given output stream.
Definition: input_output.f90:74
base_routines
This module contains all the low-level base routines e.g., all debug, control, and low-level communic...
Definition: base_routines.f90:45
coordinate_routines::coordinate_convert_from_rc
COORDINATE_CONVERT_FROM_RC performs a coordinate transformation from a rectangular cartesian coordina...
Definition: coordinate_routines.f90:119
coordinate_routines::coordinate_system_type_string
character(len=21), dimension(5), public coordinate_system_type_string
Definition: coordinate_routines.f90:106
coordinate_routines::coordinate_cylindrical_polar_type
integer(intg), parameter, public coordinate_cylindrical_polar_type
Cylindrical polar coordinate system type.
Definition: coordinate_routines.f90:71
constants::part_deriv_s1_s1
integer(intg), parameter part_deriv_s1_s1
Second partial derivative in the s1 direction i.e., d^2u/ds1ds1.
Definition: constants.f90:182
iso_varying_string::trim
Definition: iso_varying_string.f90:204
coordinate_routines::dxz
Calculates DX(:)/DZ(I) at X, where Z(I) are rectangular cartesian and X(:) are curvilinear coordinate...
Definition: coordinate_routines.f90:139
maths::matrixtranspose
Returns the transpose of a matrix A in A^T.
Definition: maths.f90:197
kinds::sp
integer, parameter sp
Single precision real kind.
Definition: kinds.f90:67
coordinate_routines::coordinates_materialsystemcalculatedxdnu2d
subroutine coordinates_materialsystemcalculatedxdnu2d(geometricInterpPointMetrics, angle, dXdNu, err, error,)
Calculates transformation between spatial CS and rotated reference orthogonal material CS in 2D space...
Definition: coordinate_routines.f90:3998
constants::part_deriv_s2_s3
integer(intg), parameter part_deriv_s2_s3
Cross derivative in the s2 and s3 direction i.e., d^2u/ds2ds3.
Definition: constants.f90:189
constants::part_deriv_s1_s3
integer(intg), parameter part_deriv_s1_s3
Cross derivative in the s1 and s3 direction i.e., d^2u/ds1ds3.
Definition: constants.f90:188
coordinate_routines::coordinate_system_destroy
subroutine, public coordinate_system_destroy(COORDINATE_SYSTEM, ERR, ERROR,)
Destroys a coordinate system.
Definition: coordinate_routines.f90:1852
coordinate_routines::coordinate_system_type_get
subroutine, public coordinate_system_type_get(COORDINATE_SYSTEM, SYSTEM_TYPE, ERR, ERROR,)
Gets the coordinate system type.
Definition: coordinate_routines.f90:1323
constants::partial_derivative_first_derivative_map
integer(intg), dimension(4) partial_derivative_first_derivative_map
PARTIAL_DERIVATIVE_FIRST_DERIVATIVE_MAP(nic) gives the partial derivative index for the first derivat...
Definition: constants.f90:254
coordinate_routines::coordinate_systems_type
Definition: coordinate_routines.f90:99
input_output::write_string_matrix_name_and_indices
integer(intg), parameter, public write_string_matrix_name_and_indices
Write the matrix name together with the matrix indices.
Definition: input_output.f90:66
coordinate_routines::coordinates_materialsystemcalculatedxdnu3d
subroutine coordinates_materialsystemcalculatedxdnu3d(geometricInterpPointMetrics, angle, dXdNu, err, error,)
Calculates transformation between spatial CS and rotated reference orthogonal material CS in 3D space...
Definition: coordinate_routines.f90:4078
input_output::writestringvalue
Write a string followed by a value to a given output stream.
Definition: input_output.f90:97
coordinate_routines::coordinate_metrics_calculate
subroutine, public coordinate_metrics_calculate(COORDINATE_SYSTEM, JACOBIAN_TYPE, METRICS, ERR, ERROR,)
Calculates the covariant metric tensor GL(i,j), the contravariant metric tensor GU(i,J), the Jacobian and derivative of the interpolated coordinate system (XI_i) with respect to the given coordinate (X_j) system (DXI_DX) at a point (X - normally a Gauss point). Old cmiss name: XGMG.
Definition: coordinate_routines.f90:648
coordinate_routines::coordinate_rectangular_cartesian_type
integer(intg), parameter, public coordinate_rectangular_cartesian_type
Rectangular Cartesian coordinate system type.
Definition: coordinate_routines.f90:70
coordinate_routines::coordinate_delta_calculate_dp
real(dp) function, dimension(size(x, 1)) coordinate_delta_calculate_dp(COORDINATE_SYSTEM, X, Y, ERR, ERROR)
Calculates the difference (or detlta) between the point X and the point Y i.e., Y-X, in the given coordinate system. 0->2Pi discontinuities with polar coordinates are accounted for.
Definition: coordinate_routines.f90:598
maths::determinant
Returns the determinant of a matrix.
Definition: maths.f90:94
coordinate_routines::coordinate_convert_to_rc
COORDINATE_CONVERT_TO_RC performs a coordinate transformation from a coordinate system identified by ...
Definition: coordinate_routines.f90:126
coordinate_routines::dzx_dp
real(dp) function, dimension(size(x, 1)) dzx_dp(COORDINATE_SYSTEM, I, X, ERR, ERROR)
Definition: coordinate_routines.f90:2306
base_routines::diagnostics1
logical, save, public diagnostics1
.TRUE. if level 1 diagnostic output is active in the current routine
Definition: base_routines.f90:153
maths::crossproduct
Calculates the vector cross product of two vectors.
Definition: maths.f90:73
coordinate_routines::d2zx
Definition: coordinate_routines.f90:145
types::field_interpolated_point_type
Contains the interpolated value (and the derivatives wrt xi) of a field at a point. Old CMISS name XG.
Definition: types.f90:1129
input_output::write_string_vector
Write a string followed by a vector to a specified output stream.
Definition: input_output.f90:290
base_routines::diagnostic_output_type
integer(intg), parameter, public diagnostic_output_type
Diagnostic output type.
Definition: base_routines.f90:67
constants::part_deriv_s1_s2
integer(intg), parameter part_deriv_s1_s2
Cross derivative in the s1 and s2 direction i.e., d^2u/ds1ds2.
Definition: constants.f90:185
constants::part_deriv_s2_s2
integer(intg), parameter part_deriv_s2_s2
Second partial derivative in the s2 direction i.e., d^2u/ds2ds2.
Definition: constants.f90:184
coordinate_routines::coordinate_system_type_set
subroutine, public coordinate_system_type_set(COORDINATE_SYSTEM, TYPE, ERR, ERROR,)
Sets/changes the type of a coordinate system.
Definition: coordinate_routines.f90:1539
coordinate_routines::coordinate_system_create_finish
subroutine, public coordinate_system_create_finish(COORDINATE_SYSTEM, ERR, ERROR,)
Finishes the creation of a new coordinate system.
Definition: coordinate_routines.f90:1810
maths::identitymatrix
Returns the identity matrix.
Definition: maths.f90:155
coordinate_routines::coordinate_delta_calculate
Calculates the difference (or delta) between two points in a coordinate system. Discontinuities for p...
Definition: coordinate_routines.f90:133
input_output::writestring
Write a string to a given output stream.
Definition: input_output.f90:80
coordinate_routines::coordinate_derivative_convert_to_rc
Definition: coordinate_routines.f90:156
constants::part_deriv_s3_s3
integer(intg), parameter part_deriv_s3_s3
Second partial derivative in the s3 direction i.e., d^2u/ds3ds3.
Definition: constants.f90:187
types::field_interpolation_parameters_type
Contains the parameters required to interpolate a field variable within an element. Old CMISS name XE.
Definition: types.f90:1141
coordinate_routines::coordinate_oblate_spheroidal_type
integer(intg), parameter, public coordinate_oblate_spheroidal_type
Oblate spheroidal coordinate system type.
Definition: coordinate_routines.f90:74
input_output::writestringvector
Write a string followed by a vector to a specified output stream.
Definition: input_output.f90:299
coordinate_routines::coordinate_system_orientation_set
subroutine, public coordinate_system_orientation_set(COORDINATE_SYSTEM, ORIENTATION, ERR, ERROR,)
Sets/changes the orientation of a coordinate system.
Definition: coordinate_routines.f90:1693
coordinate_routines::coordinate_jacobian_volume_type
integer(intg), parameter, public coordinate_jacobian_volume_type
Volume type Jacobian.
Definition: coordinate_routines.f90:94
maths::l2norm
Returns the L2 norm of a vector.
Definition: maths.f90:161
coordinate_routines::coordinate_radial_squared_interpolation_type
integer(intg), parameter, public coordinate_radial_squared_interpolation_type
r^2 radial interpolation
Definition: coordinate_routines.f90:83
base_routines::flagerror
Flags an error condition.
Definition: base_routines.f90:240
coordinate_routines::coordinate_systems_finalise
subroutine, public coordinate_systems_finalise(ERR, ERROR,)
Finalises the coordinate systems and destroys all coordinate systems.
Definition: coordinate_routines.f90:4290
maths::matrixproduct
Calculates and returns the matrix-product A*B in the matrix C.
Definition: maths.f90:173
coordinate_routines::coordinate_radial_interpolation_type
integer(intg), parameter, public coordinate_radial_interpolation_type
r radial interpolation
Definition: coordinate_routines.f90:82
constants::zero_tolerance
real(dp), parameter zero_tolerance
Definition: constants.f90:70
coordinate_routines::coordinates_radialinterpolationtypeset
subroutine, public coordinates_radialinterpolationtypeset(coordinateSystem, radialInterpolationType, err, error,)
Sets/changes the radial interpolation type of a coordinate system.
Definition: coordinate_routines.f90:1468
coordinate_routines::coordinate_system_focus_set
subroutine, public coordinate_system_focus_set(COORDINATE_SYSTEM, FOCUS, ERR, ERROR,)
Sets/changes the focus of a coordinate system.
Definition: coordinate_routines.f90:1421
kinds
This module contains all kind definitions.
Definition: kinds.f90:45
coordinate_routines::coordinate_convert_to_rc_sp
real(sp) function, dimension(size(x, 1)) coordinate_convert_to_rc_sp(COORDINATE_SYSTEM, X, ERR, ERROR)
COORDINATE_CONVERT_TO_RC_SP performs a coordinate transformation from a coordinate system identified ...
Definition: coordinate_routines.f90:522
constants::part_deriv_s1_s2_s3
integer(intg), parameter part_deriv_s1_s2_s3
Cross derivative in the s1, s2 and s3 direction i.e., d^3u/ds1ds2ds3.
Definition: constants.f90:190
coordinate_routines::coordinate_convert_from_rc_dp
real(dp) function, dimension(size(z, 1)) coordinate_convert_from_rc_dp(COORDINATE_SYSTEM, Z, ERR, ERROR)
Performs a coordinate transformation from a rectangular cartesian coordinate at the point with coordi...
Definition: coordinate_routines.f90:206
coordinate_routines::coordinate_spherical_polar_type
integer(intg), parameter, public coordinate_spherical_polar_type
Spherical polar coordinate system type.
Definition: coordinate_routines.f90:72
input_output
This module handles all formating and input and output.
Definition: input_output.f90:45