-
Notifications
You must be signed in to change notification settings - Fork 3
/
chobjmod.f90
2196 lines (2159 loc) · 79.3 KB
/
chobjmod.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
! Copyright 2021, Christian Hafner
!
! This file is part of OpenMaXwell.
!
! OpenMaXwell is free software: you can redistribute it and/or modify
! it under the terms of the GNU General Public License as published by
! the Free Software Foundation, either version 3 of the License, or
! (at your option) any later version.
! OpenMaXwell is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
! GNU General Public License for more details.
! You should have received a copy of the GNU General Public License
! along with OpenMaXwell. If not, see <http://www.gnu.org/licenses/>.
MODULE CHOBJ
! Objects
USE CHEXP
SAVE
CONTAINS
! Threads
Subroutine MTModifyObject(lCheck)
Implicit none
Include 'resource.fd'
Integer(4) idum
Logical(4), intent(in) :: lCheck
Logical(4) ldum
if(.not.lCheck) then
idum=MessageBoxQQ('Set view plane = xy plane?'C,'Modify objects'C, &
MB$YESNOCANCEL.or.MB$ICONQUESTION)
if(idum.eq.MB$IDCANCEL) return
if(idum.eq.MB$IDYES) then
viewPlane(1:3,1:3)=0.0d0 ! view plane = xy, (origin unchanged)
viewPlane(3,0)=0.0d0
viewPlane(1,1)=1.0d0
viewPlane(2,2)=1.0d0
viewPlane(3,3)=1.0d0
end if
end if
iDraBnd=0
iDraExp=0
iDomBnd=0_2
iColBnd=0_2
iConBnd=0_2
iDomExp=0_2
iColExp=0_2
iConExp=0_2
iWinAction=6_2
call MTDrawExpansion(.true.)
call MTDrawBoundary(.true.)
ldum=ModifyMenuFlagsQQ(3,1,$MenuUnChecked) !boundary
ldum=ModifyMenuFlagsQQ(3,2,$MenuUnChecked) !expansion
ldum=ModifyMenuFlagsQQ(3,3,$MenuUnChecked) !field
ldum=ModifyMenuFlagsQQ(3,4,$MenuUnChecked) !function
ldum=ModifyMenuFlagsQQ(3,5,$MenuChecked) !object
ldum=ModifyMenuFlagsQQ(3,6,$MenuUnChecked) !nothing
call OutTxt('t1','Modify objects!'C)
end Subroutine MTModifyObject
! I/O
Subroutine SaveObject(lCheck)
! save Object data in a file
Implicit none
Logical, intent(in) :: lCheck
Integer(4) iOk,ios,i,idum
if(.not.lCheck) then
call Open2write(-1,'Select Object data file to be written!','Object data file ',OBJFileName,'3DO',ios)
if(ios.gt.0) return
end if
open(1,file=OBJFileName,iostat=ios)
if(ios.ne.0) then
idum=MessageBoxQQ('Error opening file!\rCannot save data!'C,'Save Object'C, &
MB$OK.or.MB$IconExclamation)
close(1)
return
end if
call WriteStr(1,CHOBJIdent,iOK)
ich(1)=nOBJ
sch(1:8)=' Objects'
call chwrit2(1,ich,1,rch,0,sch,8,iOK)
ich(1)=iObjDra
ich(2)=0
if(lObjMat) ich(2)=1
if(lObjMatW0) ich(2)=ich(2)+2
ich(3)=0
if(lObjFlg) ich(3)=1
ich(4)=iDraObj
ich(5)=0
if(lObjHid) ich(5)=1
call chwrit2(1,ich,5,rch,0,sch,0,iOK)
rch(1)=GRChmin
rch(2)=GRChmax
call chwrit2(1,ich,0,rch,2,sch,0,iOK)
do i=1,nOBJ
ich(1)=Int4(tOBJ(i)%iTypO)
ich(2:6)=Int4(tOBJ(i)%iPar(1:5))
call chwrit2(1,ich,6,rch,0,sch,0,iOK)
rch(1)=tOBJ(i)%GrfRes
call chwrit2(1,ich,0,rch,1,sch,0,iOK)
ich(1)=Int4(tOBJ(i)%iCol)
ich(2)=Int4(tOBJ(i)%iColMin)
ich(3)=Int4(tOBJ(i)%iColMax)
call chwrit2(1,ich,3,rch,0,sch,0,iOK)
rch(1:5)=tOBJ(i)%Par(1:5)
call chwrit2(1,ich,0,rch,5,sch,0,iOK)
rch(1:3)=tOBJ(i)%Plane(1:3,0)
call chwrit2(1,ich,0,rch,3,sch,0,iOK)
rch(1:3)=tOBJ(i)%Plane(1:3,1)
call chwrit2(1,ich,0,rch,3,sch,0,iOK)
rch(1:3)=tOBJ(i)%Plane(1:3,2)
call chwrit2(1,ich,0,rch,3,sch,0,iOK)
rch(1:3)=tOBJ(i)%Plane(1:3,3)
call chwrit2(1,ich,0,rch,3,sch,0,iOK)
rch(1:3)=tOBJ(i)%O(1:3)
call chwrit2(1,ich,0,rch,3,sch,0,iOK)
rch(1:3)=tOBJ(i)%e(1:3)
call chwrit2(1,ich,0,rch,3,sch,0,iOK)
end do
sch(1:1)=' '
call chwrit2(1,ich,0,rch,0,sch,1,iOK)
EndFile(1)
close(1)
end Subroutine SaveObject
Subroutine OpenObject(lCheck)
! read Object data from a file
Implicit none
Logical, intent(in) :: lCheck
Logical ldum,lFileExist
Integer(4) iOk,ios,i,i1,idum
Character(20) text
if(.not.lgcFld) lGet3DMat=.true.
if(lAskOBJ.or.(.not.lCheck)) then
nInsObj=nOBJ
call InsertDialog(.true.)
if(kInsObj.lt.0) return
lInsertOBJ=lInsObj
end if
if(.not.lCheck) then
call Open2read(-1,'Select Object data file to be read!','Object data file ',OBJFileName,'3DO',ios)
if(ios.gt.0) return
end if
inquire(file=OBJFileName,Exist=lFileExist)
if(.not.lFileExist) return
open(1,file=OBJFileName,status='old',iostat=ios)
if(ios.ne.0) then
if(.not.lCheck) idum=MessageBoxQQ('Error opening file!\rCannot read data!'C,'Open Object'C, &
MB$OK.or.MB$IconExclamation)
close(1)
return
end if
call ReadStr(1,text,iOK)
if(CHOBJIdent(1:17).ne.text(1:17)) then
idum=MessageBoxQQ('Wrong file Type!\rContinue reading?'C,'Open Object'C, &
MB$YesNo.or.MB$IconQuestion)
if(idum.eq.MB$IDNO) then
close(1)
return
end if
end if
call chread2(1,ich,1,rch,0,iOK)
if(iOK.ne.0) then
idum=MessageBoxQQ('Error reading input file!\r(number of Objects)'C,'Open Object'C, &
MB$OK.or.MB$ICONSTOP)
close(1)
return
end if
nInsObj=ich(1)
if(lInsertOBJ) then
i=kInsObj
else
i=0
nOBJ=0
end if
call chread2(1,ich,5,rch,0,iOK)
if(iOK.ne.0) then
idum=MessageBoxQQ('Error reading input file!\r(Header data)'C,'Open Object'C, &
MB$OK.or.MB$ICONSTOP)
close(1)
return
end if
iObjDra=ich(1)
lCHGLdoubleSide=.true.
if((iObjDra.eq.3).or.(iObjDra.eq.4)) lCHGLdoubleSide=.false.
lObjMat=.false.
lObjMatW0=.false.
if((ich(2).eq.1).or.(ich(2).eq.1)) lObjMat=.true.
if(ich(2).gt.1) lObjMatW0=.true.
lObjFlg=.false.
if(ich(3).eq.1) lObjFlg=.true.
iDraObj=ich(4)
lObjHid=.false.
if(ich(5).eq.1) lObjHid=.true.
lGRCused=.not.lObjHid
call GRCclear()
call chread2(1,ich,0,rch,2,iOK)
if(iOK.ne.0) then
idum=MessageBoxQQ('Error reading input file!\r(Header data)'C,'Open Object'C, &
MB$OK.or.MB$ICONSTOP)
close(1)
return
end if
GRChmin=rch(1)
GRChmax=rch(2)
do i1=1,nInsObj
call InsertOBJ(i,1,ldum)
if(.not.ldum) then
idum=MessageBoxQQ('Cannot insert Object!'C,'Open Object'C, &
MB$OK.or.MB$ICONSTOP)
close(1)
return
end if
i=i+1
call chread2(1,ich,6,rch,0,iOK)
if(iOK.ne.0) then
idum=MessageBoxQQ('Error reading input file!\r(Object data)'C,'Open Object'C, &
MB$OK.or.MB$ICONSTOP)
close(1)
return
end if
tOBJ(i)%iTypO=Int2(ich(1))
tOBJ(i)%iPar(1:5)=Int2(ich(2:6))
call chread2(1,ich,0,rch,1,iOK)
if(iOK.ne.0) then
idum=MessageBoxQQ('Error reading input file!\r(Object data)'C,'Open Object'C, &
MB$OK.or.MB$ICONSTOP)
close(1)
return
end if
tOBJ(i)%GrfRes=rch(1)
call chread2(1,ich,3,rch,0,iOK)
if(iOK.ne.0) then
idum=MessageBoxQQ('Error reading input file!\r(Object data)'C,'Open Object'C, &
MB$OK.or.MB$ICONSTOP)
close(1)
return
end if
tOBJ(i)%iCol=Int2(ich(1))
tOBJ(i)%iColMin=Int2(ich(2))
tOBJ(i)%iColMax=Int2(ich(3))
call chread2(1,ich,0,rch,5,iOK)
if(iOK.ne.0) then
idum=MessageBoxQQ('Error reading input file!\r(Object data)'C,'Open Object'C, &
MB$OK.or.MB$ICONSTOP)
close(1)
return
end if
tOBJ(i)%Par(1:5)=rch(1:5)
call chread2(1,ich,0,rch,3,iOK)
if(iOK.ne.0) then
idum=MessageBoxQQ('Error reading input file!\r(Object data)'C,'Open Object'C, &
MB$OK.or.MB$ICONSTOP)
close(1)
return
end if
tOBJ(i)%Plane(1:3,0)=rch(1:3)
call chread2(1,ich,0,rch,3,iOK)
if(iOK.ne.0) then
idum=MessageBoxQQ('Error reading input file!\r(Object data)'C,'Open Object'C, &
MB$OK.or.MB$ICONSTOP)
close(1)
return
end if
tOBJ(i)%Plane(1:3,1)=rch(1:3)
call chread2(1,ich,0,rch,3,iOK)
if(iOK.ne.0) then
idum=MessageBoxQQ('Error reading input file!\r(Object data)'C,'Open Object'C, &
MB$OK.or.MB$ICONSTOP)
close(1)
return
end if
tOBJ(i)%Plane(1:3,2)=rch(1:3)
call chread2(1,ich,0,rch,3,iOK)
if(iOK.ne.0) then
idum=MessageBoxQQ('Error reading input file!\r(Object data)'C,'Open Object'C, &
MB$OK.or.MB$ICONSTOP)
close(1)
return
end if
tOBJ(i)%Plane(1:3,3)=rch(1:3)
call chread2(1,ich,0,rch,3,iOK)
if(iOK.ne.0) then
idum=MessageBoxQQ('Error reading input file!\r(Object data)'C,'Open Object'C, &
MB$OK.or.MB$ICONSTOP)
close(1)
return
end if
tOBJ(i)%O(1:3)=rch(1:3)
call chread2(1,ich,0,rch,3,iOK)
if(iOK.ne.0) then
idum=MessageBoxQQ('Error reading input file!\r(Object data)'C,'Open Object'C, &
MB$OK.or.MB$ICONSTOP)
close(1)
return
end if
tOBJ(i)%e(1:3)=rch(1:3)
end do
close(1)
kOBJ=nOBJ
kPar=nPar
end Subroutine OpenObject
Subroutine Open3DM(lCheck,iCM)
! read 3D matching point data from a file
Implicit none
Logical, intent(in) :: lCheck
Logical ldum,lFileExist
Integer(2) iD1,iD2
Integer(4) iOk,ios,i,i1,idum,k,n,iCM
Real(8) d
if(.not.lgcFld) lGet3DMat=.true.
if(lAskOBJ.or.(.not.lCheck)) then
nInsObj=nOBJ
call InsertDialog(.true.)
if(kInsObj.lt.0) return
lInsertOBJ=lInsObj
end if
if(.not.lCheck) then
call Open2read(-1,'Select 3D matching point data file to be read!','3DMP data file ',M3DFileName,'3DM',ios)
if(ios.gt.0) return
end if
inquire(file=M3DFileName,Exist=lFileExist)
if(.not.lFileExist) return
open(1,file=M3DFileName,status='old',iostat=ios)
if(ios.ne.0) then
if(.not.lCheck) idum=MessageBoxQQ('Error opening file!\rCannot read data!'C,'Open 3DMP'C, &
MB$OK.or.MB$IconExclamation)
close(1)
return
end if
! read number of identifiers and create a fictitious 2D boundary for each identifier
call chread2(1,ich,1,rch,0,iOK)
if(iOK.ne.0) then
idum=MessageBoxQQ('Error reading input file!\r(number of identifiers)'C,'Open 3DMP'C, &
MB$OK.or.MB$ICONSTOP)
close(1)
return
end if
n=ich(1)
i=0
do i1=1,n
call chread2(1,ich,2,rch,0,iOK)
if(iOK.ne.0) then
idum=MessageBoxQQ('Error reading input file!\r(identifiers)'C,'Open 3DMP'C, &
MB$OK.or.MB$ICONSTOP)
close(1)
return
end if
call InsertBnd(i,1,ldum)
if(.not.ldum) then
idum=MessageBoxQQ('Cannot insert boundary!'C,'Open 3DMP'C, &
MB$OK.or.MB$ICONSTOP)
close(1)
return
end if
i=i+1
tBnd(i)%nEdge=0
tBnd(i)%iRDom=Int2(ich(1))
tBnd(i)%iLDom=Int2(ich(2))
tBnd(i)%iCond=0_2 ! default continuity condition
call GetBndBC(i)
tBnd(i)%iConn=0_2
tBnd(i)%iTypB=2_2
tBnd(i)%iCol=Int2(iCM)
tBnd(i)%nMatPts=1
tBnd(i)%Weight=1.0d0
tBnd(i)%val(1)=0.0d0
tBnd(i)%val(2)=0.0d0
tBnd(i)%Weight2=tBnd(i)%Weight
tBnd(i)%nSpline=0
tBnd(i)%fAmpl=0.0d0
tBnd(i)%Formula='0'C
call InsertBndEdg(i,0,1,ldum)
if(.not.ldum) then
idum=MessageBoxQQ('Cannot insert corner!'C,'Open 3DMP'C, &
MB$OK.or.MB$ICONSTOP)
close(1)
return
end if
k=tBnd(i)%iEdgeOffset+1
tBndEdg(k)%x=1.0d10 ! dummy large circle
tBndEdg(k)%y=1.0d10
tBndEdg(k)%r=1.0d9
tBndEdg(k)%xa=tBndEdg(k)%x+tBndEdg(k)%r
tBndEdg(k)%ya=tBndEdg(k)%y
tBndEdg(k)%xb=tBndEdg(k)%x+tBndEdg(k)%r
tBndEdg(k)%yb=tBndEdg(k)%y
tBndEdg(k)%xo=tBndEdg(k)%x
tBndEdg(k)%yo=tBndEdg(k)%y
end do
! now read 3D matching points and convert to rectangles
call chread2(1,ich,1,rch,0,iOK)
if(iOK.ne.0) then
idum=MessageBoxQQ('Error reading input file!\r(Number of 3D matching points)'C,'Open 3DMP'C, &
MB$OK.or.MB$ICONSTOP)
close(1)
return
end if
n=ich(1)
if(lInsertOBJ) then
i=kInsObj
else
i=0
nOBJ=0
end if
do i1=1,n
call InsertOBJ(i,1,ldum)
if(.not.ldum) then
idum=MessageBoxQQ('Cannot insert Object!'C,'Open 3DMP'C, &
MB$OK.or.MB$ICONSTOP)
close(1)
return
end if
i=i+1
tOBJ(i)%iTypO=4_2 ! rectangle
call chread2(1,ich,1,rch,6,iOK)
if(iOK.ne.0) then
idum=MessageBoxQQ('Error reading input file!\r(3D matching point data)'C,'Open 3DMP'C, &
MB$OK.or.MB$ICONSTOP)
close(1)
return
end if
tOBJ(i)%iPar(1)=Int2(ich(1))
iD1=tBnd(ich(1))%iRDom
iD2=tBnd(ich(1))%iLDom
tOBJ(i)%iPar(2:5)=0_2
tOBJ(i)%iCol=Int2(iCM)
tOBJ(i)%iColMin=30_2
tOBJ(i)%iColMax=90_2
tOBJ(i)%Plane(1:3,0)=rch(1:3)
tOBJ(i)%Plane(1:3,1)=rch(4:6) ! this is wrong direction!
tOBJ(i)%Plane(1:3,2)=rch(4:6) ! this is wrong direction!
tOBJ(i)%Plane(1:3,3)=rch(4:6) ! this is correct direction!
call Ortho3DSpace3(tOBJ(i)%Plane(1:3,3),tOBJ(i)%Plane(1:3,1),tOBJ(i)%Plane(1:3,2)) ! correct wrong directions (perpendicular)
d=dsqrt(r3Vec_Length(rch(4:6)))
tOBJ(i)%GrfRes=d/10.0d0
tOBJ(i)%Plane(1:3,0)=tOBJ(i)%Plane(1:3,0)-0.5d0*d*tOBJ(i)%Plane(1:3,1)-0.5d0*d*tOBJ(i)%Plane(1:3,2) ! origin of rectangle -> 3DMP in center
tOBJ(i)%Par(1:5)=d ! dx,dy,mPsideLength,...
tOBJ(i)%Par(4)=1.0d0 ! multipole density
tOBJ(i)%O(1:3)=0.0d0
tOBJ(i)%e(1)=1.0d0
tOBJ(i)%e(2:3)=0.0d0
end do
close(1)
kOBJ=nOBJ
kPar=nPar
end Subroutine Open3DM
Subroutine Open3DD(lCheck,nOb,iTyD,iCD)
! read 3D dipole data from a file
Implicit none
Logical, intent(in) :: lCheck
Logical ldum,lFileExist
Integer(4) iOk,ios,i,i1,idum,n,nOb,iTyD,iCD
if(lAskExp.or.(.not.lCheck)) then
nInsObj=nExp
call InsertDialog(.true.)
if(kInsObj.lt.0) return
lInsertExp=lInsObj
end if
if(.not.lCheck) then
call Open2read(-1,'Select 3D dipole data file to be read!','Expansion data file ',D3DFileName,'EXP',ios)
if(ios.gt.0) return
end if
inquire(file=D3DFileName,Exist=lFileExist)
if(.not.lFileExist) return
open(1,file=D3DFileName,status='old',iostat=ios)
if(ios.ne.0) then
if(.not.lCheck) idum=MessageBoxQQ('Error opening file!\rCannot read data!'C,'Open 3D dipoles'C, &
MB$OK.or.MB$IconExclamation)
close(1)
return
end if
! read number of 3D dipoles
call chread2(1,ich,1,rch,0,iOK)
if(iOK.ne.0) then
idum=MessageBoxQQ('Error reading input file!\r(number of 3D dipoles)'C,'Open 3DD'C, &
MB$OK.or.MB$ICONSTOP)
close(1)
return
end if
n=ich(1)
if(lInsertExp) then
i=kInsObj
else
i=0
nExp=0
nPar=0
end if
do i1=1,n
call chread2(1,ich,1,rch,3,iOK)
if(iOK.ne.0) then
idum=MessageBoxQQ('Error reading input file!\r(3D dipole data)'C,'Open 3DD'C, &
MB$OK.or.MB$ICONSTOP)
close(1)
return
end if
call InsertExp(i,1,ldum)
if(.not.ldum) then
idum=MessageBoxQQ('Cannot insert expansion!'C,'Open 3DD'C, &
MB$OK.or.MB$ICONSTOP)
close(1)
return
end if
i=i+1
tExp(i)%iDom=Int2(ich(1))
tExp(i)%iTypE=6_2 ! dipole is a 3D multipole
tExp(i)%iHE=Min(2_2,Max(0_2,Int2(iTyD))) ! only electric dipoles in Uli Koch's *.3dd file -> iTyD=0
tExp(i)%iConn=0_2
tExp(i)%iCol=Int2(iCD)
tExp(i)%iObj=Int2(nOb)
tExp(i)%iE(1:5)=1
tExp(i)%iE(3)=0
tExp(i)%iE(6)=0
tExp(i)%rE(1:5)=0.0d0
tExp(i)%depend=1.0d100
tExp(i)%gc=(0.0d0,0.0d0)
tExp(i)%xo=rch(1)
tExp(i)%yo=rch(2)
tExp(i)%Plane(1:3,0)=rch(1:3)
tExp(i)%Plane(1:3,1:3)=0.0d0
tExp(i)%Plane(1,1)=1.0d0
tExp(i)%Plane(2,2)=1.0d0
tExp(i)%Plane(3,3)=1.0d0
tExp(i)%O(1:3)=0.0d0
tExp(i)%e(1:3)=0.0d0
tExp(i)%nPar=0
if(tExp(i)%iHE.ne.2_2) then
call InsertPar(i,0,3,ldum) ! electric and magnetic dipoles have 3 parameters
else
call InsertPar(i,0,6,ldum) ! electric + magnetic dipoles have 6 parameters
end if
if(.not.ldum) then
idum=MessageBoxQQ('Cannot insert parameter!'C,'Open expansion'C, &
MB$OK.or.MB$ICONSTOP)
close(1)
return
end if
end do
close(1)
call CorrExpPar(1000_4)
kExp=nExp
kPar=nPar
end Subroutine Open3DD
! Object operations
Subroutine gen3DMatPts(k1,k2,lcount,iAlsoW,nPts,Points,iC,iEQ,iOC)
! generate 3D matching points for the 3D objects k1...k2 -> Points
! if iC is present and iC(1)=0, save the color number for the points in iC
! if iC is present and iC(1)=1, save the number of the original 2D matching point in iC
! if iC is present and iC(1)=2, save the domain numbers of the original 2D matching point in iC
! if iC is present and iC(1)=3, save the boundary number of the original 2D matching point in iC
! if iEQ is present: save number of equations in iEQ
! if iOC is present: save number of 3D object in iOC
! nPts stores the numbers of matching points for each object
! iAlsoW=0: ignore boundaries with weight 0;
! 1: use boundaries with weight 0 if domain numbers are different
! 2: use all boundaries
Implicit none
Integer(4) k1,k2,kO,k,i1,i2,i,nPts(*),nP,mP,idum,nEq2D,na,nha,nBndPtI1,nBndPtI2
Integer(2) iDl,iDR,iCl(1)
Integer(2), Optional:: iC(*),iEQ(*),iOC(*)
Real(8) Points(3,0:3,*),P(3),vt(3),s,ds,r0,PA(2),PB(2),PC(2),a,b,c,F,ha,h,dh,xa,x,dx
Logical lcount,l2DBnd,l2DDom,lsav,l2DNr,lAlsoW
Integer(4), intent(in) :: iAlsoW
if(((iBound.ne.2_2).and.(iBound.ne.4_2))) call cBndGetABO() ! get c-poly, splines, match.pts
nBndEq3D=0
nBndPt3D=0
lsav=.false.
l2DBnd=.false.
l2DDom=.false.
l2DNr=.false.
if(Present(iC)) then
lsav=.true.
l2DBnd=.true.
l2DDom=.true.
if(iC(1).eq.0_2) l2DBnd=.false.
if(iC(1).eq.1_2) l2DDom=.false.
if(iC(1).eq.3_2) l2DNr=.true.
end if
i1=0
do kO=k1,k2
mP=0
if((tOBJ(kO)%iTypO.eq.3).or.(tOBJ(kO)%iTypO.eq.4)) then ! Triangle or Rectangle
if(tBnd(max(1,tOBJ(kO)%iPar(1)))%nMat.lt.1) Cycle ! no matching points for this object
nEq2D=0
iDR=tBnd(tOBJ(kO)%iPar(1))%iLDom
iDL=tBnd(tOBJ(kO)%iPar(1))%iRDom
lAlsoW=.true.
if((iAlsoW.lt.1).or.((iAlsoW.eq.1).and.(iDR.eq.IDL))) lAlsoW=.false.
if((iDR.lt.-30000_2).or.((tBnd(tOBJ(kO)%iPar(1))%weight.lt.pSmall)).and.(.not.lAlsoW)) Cycle
nEq2D=nEq2D+tBnd(tOBJ(kO)%iPar(1))%nBC
! if(iDL.ne.0_2) nEq2D=nEq2D+3
! if(iDR.ne.0_2) nEq2D=nEq2D+3
if(tOBJ(kO)%iTypO.eq.3) then ! Triangle
PB(1:2)=tOBJ(kO)%O(1:2)
PC(1)=tOBJ(kO)%O(3)
PC(2)=tOBJ(kO)%e(1)
PA(1:2)=tOBJ(kO)%e(2:3)
call getTriangle2DData(PA,PB,PC,a,b,c,s,F,ha)
nP=0
nha=min(100_4,max(1_4,nint(ha/tOBJ(kO)%Par(1),4)))
dh=ha/Dble(nha)
h=-0.5d0*dh
do i=1,nha
h=h+dh
na=min(100_4,max(1_4,nint(a*(1.0d0-h/ha)/tOBJ(kO)%Par(1),4)))
nP=nP+na
end do
else ! Rectangle
na= min(100_4,max(1_4,nint(tOBJ(kO)%Par(1)/tOBJ(kO)%Par(3),4)))
nha=min(100_4,max(1_4,nint(tOBJ(kO)%Par(2)/tOBJ(kO)%Par(3),4)))
dx=tOBJ(kO)%Par(1)/Dble(na)
dh=tOBJ(kO)%Par(2)/Dble(nha)
nP=na*nha
end if
mP=mP+nP
nBndEq3D=nBndEq3D+nEq2D*nP
nBndPt3D=nBndPt3D+nP
if(.not.lcount) then
i2=i1+nP
i1=i1+1
if(iDL.ne.-257_2) then
if(iObjDra.gt.1) then ! use domain colors
iC(i1)=iPackDom(nDom,iDL,iDR)
if(i2.gt.i1) iC(i1+1:i2)=iC(i1)
else
iC(i1:i2)=-30001_2
end if
end if
i=tBnd(tOBJ(kO)%iPar(1))%iMatOffset+1
if(lsav.and.l2DBnd.and.(.not.l2DDom)) iC(i1:i2)=Int2(i)
if(l2DNr) iC(i1:i2)=Int2(tOBJ(kO)%iPar(1))
if(Present(iEQ)) iEQ(i1:i2)=nEq2D
if(Present(iOC)) iOC(i1:i2)=kO
Points(1:2,3,i1)=0.0d0
Points(3,3,i1)=dh**2
if(tOBJ(kO)%iTypO.eq.3) then ! Triangle
Points(1:2,1,i1)=PC(1:2)-PB(1:2)
call Unit2DV(Points(1:2,1,i1))
Points(1:2,1,i1)=tOBJ(kO)%Par(1)*Points(1:2,1,i1)
Points(3,1:2,i1)=0.0d0
Points(1,2,i1)=-Points(2,1,i1)
Points(2,2,i1)=Points(1,1,i1)
call rvLoc2Glob(Points(1:3,1,i1),tObj(kO)%Plane,Points(1:3,1,i1))
call rvLoc2Glob(Points(1:3,2,i1),tObj(kO)%Plane,Points(1:3,2,i1))
call rvLoc2Glob(Points(1:3,3,i1),tObj(kO)%Plane,Points(1:3,3,i1))
do i=i1+1,i2
Points(1:3,0:3,i)=Points(1:3,0:3,i1)
end do
i1=i1-1
h=-0.5d0*dh
do i=1,nha
h=h+dh
xa=a*(1.0d0-h/ha)
na=min(100_4,max(1_4,nint(xa/tOBJ(kO)%Par(1),4)))
dx=xa/Dble(na)
x=-0.5d0*dx
xa=h*dsqrt(dabs((c/ha)**2-1.0d0))
do k=1,na
x=x+dx
i1=i1+1
Points(1,0,i1)=xa+x
Points(2,0,i1)=h
Points(3,0,i1)=0.0d0
call vLoc2Glob(Points(1:3,0,i1),tObj(kO)%Plane,Points(1:3,0,i1))
end do
end do
else ! Rectangle
Points(1:3,1:2,i1)=0.0d0
Points(1,1,i1)=dx
Points(2,2,i1)=dh
call rvLoc2Glob(Points(1:3,1,i1),tObj(kO)%Plane,Points(1:3,1,i1))
call rvLoc2Glob(Points(1:3,2,i1),tObj(kO)%Plane,Points(1:3,2,i1))
call rvLoc2Glob(Points(1:3,3,i1),tObj(kO)%Plane,Points(1:3,3,i1))
do i=i1+1,i2
Points(1:3,0:3,i)=Points(1:3,0:3,i1)
end do
i1=i1-1
h=-0.5d0*dh
do i=1,nha
h=h+dh
x=-0.5d0*dx
do k=1,na
x=x+dx
i1=i1+1
Points(1,0,i1)=x
Points(2,0,i1)=h
Points(3,0,i1)=0.0d0
call vLoc2Glob(Points(1:3,0,i1),tObj(kO)%Plane,Points(1:3,0,i1))
end do
end do
end if
end if
nPts(kO-k1+1)=mP
Cycle
end if
if((tOBJ(kO)%iPar(1).gt.tOBJ(kO)%iPar(2)).or.(tOBJ(kO)%iPar(2).lt.1)) Cycle
do k=max(1,tOBJ(kO)%iPar(1)),tOBJ(kO)%iPar(2)
nBndPtI1=tBnd(k)%iMatOffset+1
nBndPtI2=tBnd(k)%iMatOffset+tBnd(k)%nMat
do i=nBndPtI1,nBndPtI2
s=sBndPt(i)
call GetBndPt(k,s,P(1),P(2),vt(1),vt(2),iDL,iDR,idum)
nEq2D=0
lAlsoW=.true.
if((iAlsoW.lt.1).or.((iAlsoW.eq.1).and.(iDR.eq.IDL))) lAlsoW=.false.
if((iDL.lt.-30000_2).or.((tBnd(idum)%weight.lt.pSmall)).and.(.not.lAlsoW)) Cycle
nEq2D=nEq2D+tBnd(k)%nBC
! if(iDL.ne.0_2) nEq2D=nEq2D+3
! if(iDR.ne.0_2) nEq2D=nEq2D+3
P(3)=0.0d0
vt(3)=0.0d0
if(i.gt.nBndPtI1) then
if(i.lt.nBndPtI2) then
ds=0.5d0*(sBndPt(i+1)-sBndPt(i-1))
else
ds=sBndPt(i)-sBndPt(i-1)
end if
else if(nBndPtI2.gt.nBndPtI1) then
ds=sBndPt(i+1)-sBndPt(i)
else
ds=0.5d0*tBnd(k)%sLength
end if
vt=ds*vt
Select Case(tOBJ(kO)%iTypO)
Case(0) ! Torus
s=P(1)*Pi*abs(tOBJ(kO)%Par(2))/180.0d0
Case(1) ! Cylinder
s=abs(tOBJ(kO)%Par(2))
Case(2) ! Spiral
r0=Dist3DPtAxis(P,tOBJ(kO)%O(1:3),tOBJ(kO)%e(1:3))
s=SpiralLength(r0,tOBJ(kO)%Par(1),tOBJ(kO)%Par(2),tOBJ(kO)%Par(3))
Case(5) ! Cone
s=abs(tOBJ(kO)%Par(2))
end Select
if(tOBJ(kO)%Par(4).le.0.0d0) then ! fixed number of points
nP=max(1_4,min(10000_4,nint(-tOBJ(kO)%Par(4),4)))
else
nP=max(1_4,min(10000_4,nint(tOBJ(kO)%Par(4)*abs(s.div.ds),4)))
end if
mP=mP+nP
nBndEq3D=nBndEq3D+nEq2D*nP
nBndPt3D=nBndPt3D+nP
if(.not.lcount) then
i2=i1+nP
i1=i1+1
Select Case(tOBJ(kO)%iTypO)
Case(0) ! Torus
if(lsav.and.((.not.l2DBnd).or.l2DDom)) then
call Pt2Torus(P,tOBJ(kO)%Par(1),tOBJ(kO)%Par(2),-nP,tOBJ(kO)%Plane,Points(1:3,0,i1:i2), &
& iObjDra,nDom,iDL,iDR,iC(i1:i2),Points(1:3,1,i1:i2),Points(1:3,2,i1:i2),Points(1:3,3,i1:i2),vt)
else
call Pt2Torus(P,tOBJ(kO)%Par(1),tOBJ(kO)%Par(2),-nP,tOBJ(kO)%Plane,Points(1:3,0,i1:i2), &
& iObjDra,nDom,-257_2,iDR,iCl,Points(1:3,1,i1:i2),Points(1:3,2,i1:i2),Points(1:3,3,i1:i2),vt)
end if
Case(1) ! Cylinder
if(lsav.and.((.not.l2DBnd).or.l2DDom)) then
call Pt2Cylinder(P,tOBJ(kO)%Par(1),tOBJ(kO)%Par(2),-nP,tOBJ(kO)%Plane,Points(1:3,0,i1:i2), &
& iObjDra,nDom,iDL,iDR,iC(i1:i2),Points(1:3,1,i1:i2),Points(1:3,2,i1:i2),Points(1:3,3,i1:i2),vt)
else
call Pt2Cylinder(P,tOBJ(kO)%Par(1),tOBJ(kO)%Par(2),-nP,tOBJ(kO)%Plane,Points(1:3,0,i1:i2), &
& iObjDra,nDom,-257_2,iDR,iCl,Points(1:3,1,i1:i2),Points(1:3,2,i1:i2),Points(1:3,3,i1:i2),vt)
end if
Case(2) ! Spiral
if(lsav.and.((.not.l2DBnd).or.l2DDom)) then
call Pt2Spiral(P,tOBJ(kO)%O,tOBJ(kO)%e,tOBJ(kO)%Par(1),tOBJ(kO)%Par(2),tOBJ(kO)%Par(3), &
& -nP,tOBJ(kO)%Plane,Points(1:3,0,i1:i2), &
& iObjDra,nDom,iDL,iDR,iC(i1:i2),Points(1:3,1,i1:i2),Points(1:3,2,i1:i2),Points(1:3,3,i1:i2),vt)
else
call Pt2Spiral(P,tOBJ(kO)%O,tOBJ(kO)%e,tOBJ(kO)%Par(1),tOBJ(kO)%Par(2),tOBJ(kO)%Par(3), &
& -nP,tOBJ(kO)%Plane,Points(1:3,0,i1:i2), &
& iObjDra,nDom,-257_2,iDR,iCl,Points(1:3,1,i1:i2),Points(1:3,2,i1:i2),Points(1:3,3,i1:i2),vt)
end if
Case(5) ! Cone
if(lsav.and.((.not.l2DBnd).or.l2DDom)) then
call Pt2Cone(P,tOBJ(kO)%O,tOBJ(kO)%Par(1),tOBJ(kO)%Par(2),-nP,tOBJ(kO)%Plane,Points(1:3,0,i1:i2), &
& iObjDra,nDom,iDL,iDR,iC(i1:i2),Points(1:3,1,i1:i2),Points(1:3,2,i1:i2),Points(1:3,3,i1:i2),vt)
else
call Pt2Cone(P,tOBJ(kO)%O,tOBJ(kO)%Par(1),tOBJ(kO)%Par(2),-nP,tOBJ(kO)%Plane,Points(1:3,0,i1:i2), &
& iObjDra,nDom,-257_2,iDR,iCl,Points(1:3,1,i1:i2),Points(1:3,2,i1:i2),Points(1:3,3,i1:i2),vt)
end if
end Select
if(lsav.and.l2DBnd.and.(.not.l2DDom)) iC(i1:i2)=Int2(i)
if(l2DNr) iC(i1:i2)=Int2(k)
if(Present(iEQ)) iEQ(i1:i2)=nEq2D
if(Present(iOC)) iOC(i1:i2)=kO
i1=i2
end if
end do
end do
nPts(kO-k1+1)=mP
end do
end Subroutine gen3DMatPts
Subroutine get3DMatPts(k1,k2,iAlsoW,iWhat0,lInhibit)
! count 3D matching points, allocate memory, compute
! iWhat=0: iBndPt3D contains color number
! iWhat=1: iBndPt3D contains number of corresponding 2D matching point
! iWhat=2: iBndPt3D contains domain numbers of corresponding 2D matching point
! iWhat=3: iBndPt3D contains boundary number of corresponding 2D matching point
! iWhat<0: use same as before
! iAlsoW=0: ignore boundaries with weight 0;
! 1: use boundaries with weight 0 if domain numbers are different
! 2: use all boundaries
Implicit none
Real(8) Pdum(3)
Integer(4) k1,k2,kO,ier,idum,kPt
Integer(2) iWhat0,iWhat
Integer(4), intent(in) :: iAlsoW
Logical, intent(in) :: lInhibit
iWhat=iWhat0
if(iWhat.lt.0) iWhat=iWhatiBndPt3D
if((iWhat.eq.iWhatiBndPt3D).and.(.not.lGet3DMat)) return
! count 3D matching points
tOBJ(1:nObj)%nMat=0
tOBJ(1:nObj)%iMatOffset=0
call gen3DMatPts(k1,k2,.true.,iAlsoW,tOBJ(k1:k2)%nMat,Pdum)
! allocate memory
if(.not.lStopThread) then
if((.not.Allocated(eBndPt3D)).or.(nBndPt3D.ne.mBndPt3DAlloc)) then
if(Allocated(eBndPt3D)) DeAllocate(eBndPt3D)
if(Allocated(fBndPt3D)) DeAllocate(fBndPt3D)
if(Allocated(BndPt3D)) Deallocate(BndPt3D)
if(Allocated(iBndPt3D)) Deallocate(iBndPt3D)
if(Allocated(iObjBndPt3D)) Deallocate(iObjBndPt3D)
if(Allocated(iEquBndPt3D)) Deallocate(iEquBndPt3D)
idum=max(1_4,nBndPt3D)
Allocate(eBndPt3D(idum),fBndPt3D(idum),BndPt3D(3,0:3,idum),iBndPt3D(idum),iObjBndPt3D(idum),&
& iEquBndPt3D(idum),stat=ier)
if(ier.ne.0) then
idum=MessageBoxQQ('Memory alloction failed!'C,'Setup MMP matrix'C, &
MB$OK.or.MB$ICONSTOP)
mBndPt3DAlloc=-1
return
end if
mBndPt3DAlloc=idum
eBndPt3D=0.0d0
fBndPt3D=1.0d0
iWhatiBndPt3D=-1_2
end if
end if
! compute 3D matching points, count max. number of equations, set offset info
if(.not.lStopThread) then
iBndPt3D(1)=iWhat
call gen3DMatPts(k1,k2,.false.,iAlsoW,tOBJ(k1:k2)%nMat,BndPt3D,iBndPt3D,iEquBndPt3D,iObjBndPt3D)
tOBJ(k1)%iMatOffset=0
do kO=k1+1,k2
tOBJ(kO)%iMatOffset=tOBJ(kO-1)%iMatOffset+tOBJ(kO-1)%nMat
end do
end if
iWhatiBndPt3D=iWhat
lGet3DMat=.false.
! inhibit 3D matching points
if(lInhibit) then
call InhibitMatPts()
nInhBndPt3D=0
nInhEqu3D=0
do kPt=1,nBndPt3D
if(iObjBndPt3D(kPt).lt.0) then
nInhBndPt3D=nInhBndPt3D+1
nInhEqu3D=nInhEqu3D+iEquBndPt3D(kPt)
end if
end do
end if
end Subroutine get3DMatPts
Subroutine GetExpsNearMatPts(k1,k2)
! count number of expansions associated with 3D matching point and store in iBndPt3D
! for objects k1...k2
Real(8) f,fminL,fminR,f0
Integer(4) k1,k2,k,kPt,kEx
Integer(2) iL,iR,iD
if(lgcFld) return
f0=1.4d0
do k=k1,k2 ! loop over 3D objects
do kPt=tOBJ(k)%iMatOffset+1,tOBJ(k)%iMatOffset+tOBJ(k)%nMat ! loop over matching points
if(iObjBndPt3D(kPt).lt.0) Cycle ! ignore dummy object (for integrals)
iL=0
iR=0
fminL=1.0d300
fminR=1.0d300
do kEx=1,nExp
if(LSkipExp(kEx)) Cycle
if(igetiHE2(iHEGlobal,kEx).lt.0_2) Cycle
if(tExp(kEx)%iObj.lt.0) Cycle
iD=tExp(kEx)%iDom
if((iD.gt.0).and.(iD.ne.tBnd(iBndPt3D(kPt))%iRDom).and.(iD.ne.tBnd(iBndPt3D(kPt))%iLDom)) Cycle
Select Case(tExp(kEx)%iTypE)
Case(1_2) ! 2D multipole for 3D
Case(6_2) ! 3D multipole
f=r3Vec_Length(tExp(kEx)%Plane(1:3,0)-BndPt3D(1:3,0,kPt)).div.tExp(kEx)%dmin
if(iD.eq.tBnd(iBndPt3D(kPt))%iLDom) then
if(f.lt.f0) iL=min(215_2,iL+1_2)
if(f.lt.fminL) fminL=f
else if(iD.eq.tBnd(iBndPt3D(kPt))%iRDom) then
if(f.lt.f0) iR=min(215_2,iR+1_2)
if(f.lt.fminR) fminR=f
end if
Case(8_2) ! 3D ring multipole
Case(9_2) ! 3D line multipole
Case(10_2) ! 3D spiral multipole
Case(12_2) ! 2D multilayer
Case(13_2) ! 3D multilayer
end Select
end do
if(iL.gt.3_2) then
iL=5_2
else if(iL.eq.3_2) then
iL=2_2
else if(iL.eq.2_2) then
iL=3_2
else if(iL.eq.1_2) then
iL=4_2
else
iL=210_2-nint(Max(Min((20.0d0*fminL/f0),193.0d0),20.0d0),2) ! color numbers 190-17
end if
if(iR.gt.3_2) then
iR=5_2
else if(iR.eq.3_2) then
iR=2_2
else if(iR.eq.2_2) then
iR=3_2
else if(iR.eq.1_2) then
iR=4_2
else
iR=210_2-nint(Max(Min((20.0d0*fminR/f0),193.0d0),20.0d0),2)
end if
iBndPt3D(kPt)=iPackObjC(iL,iR)
end do
end do
end Subroutine GetExpsNearMatPts
Subroutine DistPtObj(kO,lO,r,lNoPer,dmin0,rNmin0,iDom0,val0,lInh,lRSidemin0,iDom2)
! find the boundary point rNmin0 with the shortest distance dmin0 from the given
! point r. Check objects kO...lO (all objects except object lO for kO=0,...)
! return domain number iDom0, and boundary values val0 in addition to dmin0,rNmin0
Implicit none
Real(8) r(3),dmin,dmi,rNmin(3),rNmi(3),dmin0,rNmin0(3),val(2),val0(2),rP0(3),rP(3),fv,d,f0,r0,x0,f1,f2,fNmin,a,b
Integer(4) k1,k2,k,kPtmin,kPt,ixzSymm0,iyzSymm0
Integer(2), Intent(in):: kO,lO
Integer(2) iDom,iDom0,iTO,iDom1
Integer(2), Optional:: iDom2
Logical, intent(in) :: lNoPer,lInh
Logical lLoc,lRSidemin
Logical, Optional:: lRSidemin0
if(lgcFld) then
call DistPtBnd(kO,lO,r(1),r(2),.true.,lNoPer,dmin0,rNmin0(1),rNmin0(2),lRSidemin,iDom0,val0,iDom1)
if(Present(lRSidemin0)) lRSidemin0=lRSidemin
if(Present(iDom2)) iDom2=iDom1
rNmin0(3)=r(3)
return
end if
if((kO.ge.1).and.(kO.le.nObj)) then
if((lO.ge.1).and.(lO.le.nObj)) then
k1=min(kO,lO)
k2=max(kO,lO)
else
k1=kO
k2=k1
end if
else if((kO.lt.0).and.(-kO.lt.nObj)) then
k1=1
k2=-kO
else
k1=1
k2=nObj
end if
rP0=r
fv=1.0d0
if(iyzSymm.ne.0) rP0(1)=dabs(rP0(1))
if(ixzSymm.ne.0) rP0(2)=dabs(rP0(2))
if(ixySymm.ne.0) rP0(3)=dabs(rP0(3))
if((iyzSymm.eq.2).and.(dabs(r(1)-rP0(1)).gt.pSmall)) fv=-fv
if((ixzSymm.eq.2).and.(dabs(r(2)-rP0(2)).gt.pSmall)) fv=-fv
if((ixySymm.eq.2).and.(dabs(r(3)-rP0(3)).gt.pSmall)) fv=-fv
call MoveToOrigCell(rP0)
dmin0=pBig
rNmin0(1:3)=0.0d0
iDom0=-9_2
iDom1=-9_2
val0(1)=0.0d0
val0(2)=0.0d0
do k=k1,k2
if((kO.eq.0).and.(lO.eq.k)) Cycle
if((tOBJ(k)%iTypO.ne.3).and.(tOBJ(k)%iTypO.ne.4).and.((tOBJ(k)%iPar(1).gt.tOBJ(k)%iPar(2)) &
& .or.(tOBJ(k)%iPar(2).lt.1))) Cycle
iTO=tOBJ(k)%iTypO
if(lInh.and.(tOBJ(k)%nInhibited.gt.0)) then
if(tOBJ(k)%nInhibited.ge.tOBJ(k)%nMat) Cycle
iTO=-1
end if
Select Case(iTO)